patch RATTLE to zero velocities on 4-site water extra sites

jayponder · jayponder · commit d88d4d930f7b · 2025-02-18T21:47:52.000-05:00
diff --git a/example/mm3.xyz b/example/mm3.xyz
@@ -4,4 +4,4 @@
      3  H      0.105212    0.000000   -0.940767    21     1
      4  O     -2.940698    0.000000   -0.082189     6     5     6
      5  H     -3.491105   -0.749216    0.097557    21     4
-     6  H     -3.491105    0.749217    0.097557    21     4
+     6  H     -3.491105    0.749216    0.097557    21     4
diff --git a/example/tip4pbox.key b/example/tip4pbox.key
@@ -1,8 +1,34 @@
-parameters     ../params/oplsaa08
-#verbose
+parameters              NONE
+rattle                  WATER
 
-x-axis                       31.1
+vdwtype                 LENNARD-JONES
+radiusrule              GEOMETRIC
+radiustype              SIGMA
+radiussize              DIAMETER
+epsilonrule             GEOMETRIC
+dielectric              1.0000
+
+x-axis                  31.125
 ewald
-ewald-cutoff                  9.0
-vdw-cutoff                   12.0
-rattle                      WATER
+ewald-cutoff            9.0
+vdw-cutoff              12.0
+vdw-correlation
+
+## TIP4P-Ew Parameters (J Chem Phys, 120, 9665-9678, 2004)
+
+atom             1    O       "TIP4P-Ew O"     8     15.995     3
+atom             2    H       "TIP4P-Ew H"     1      1.008     1
+atom             3    M       "TIP4P-Ew M"     0      0.000     1
+
+bond             1    2            1000.0      0.9572
+bond             1    3            1000.0      0.1250
+angle            2    1    2        100.0      104.52
+angle            2    1    3        100.0       52.26
+
+vdw              1                3.16435      0.162750
+vdw              2                  0.000      0.0000
+vdw              3                  0.000      0.0000
+
+charge           1                  0.000
+charge           2                0.52422
+charge           3               -1.04844
diff --git a/example/tip4pbox.xyz b/example/tip4pbox.xyz
diff --git a/source/control.f b/source/control.f
@@ -45,7 +45,7 @@ subroutine control
 c     check for control parameters on the command line
 c
       exist = .false.
-      do i = 1, narg-1
+      do i = 1, narg
          string = arg(i)
          call upcase (string)
          if (string(1:2) .eq. '-D') then
diff --git a/source/mdinit.f b/source/mdinit.f
@@ -185,7 +185,7 @@ subroutine mdinit (dt)
             if (atomic(i) .le. 0)  ndummy = ndummy + 1
             if (use(i) .and. mass(i).le.0.0d0) then
                mass(i) = 1.0d0
-               totmass = totmass + 1.0d0
+               if (atomic(i) .gt. 0)  totmass = totmass + 1.0d0
                if (verbose) then
                   write (iout,30)  i
    30             format (/,' MDINIT  --  Warning, Mass of Atom',
diff --git a/source/rattle.f b/source/rattle.f
@@ -150,7 +150,7 @@ subroutine rattle (dt,xold,yold,zold)
      &              ' Iterations')
       end if
 c
-c     set four-site water extra centers from H-O-H coordinates
+c     set position and velocity of four-site water extra centers
 c
       do i = 1, nwat4
          ia = iwat4(1,i)
@@ -162,6 +162,9 @@ subroutine rattle (dt,xold,yold,zold)
          x(ia) = oterm*x(ib) + hterm*(x(ic)+x(id))
          y(ia) = oterm*y(ib) + hterm*(y(ic)+y(id))
          z(ia) = oterm*z(ib) + hterm*(z(ic)+z(id))
+         v(1,ia) = 0.0d0
+         v(2,ia) = 0.0d0
+         v(3,ia) = 0.0d0
       end do
 c
 c     apply group position and velocity constraints via exact reset
@@ -522,7 +525,7 @@ subroutine shake (xold,yold,zold)
      &              ' Iterations')
       end if
 c
-c     set four-site water extra centers from H-O-H coordinates
+c     set the position of any four-site water extra centers
 c
       do i = 1, nwat4
          ia = iwat4(1,i)
@@ -719,7 +722,7 @@ subroutine water4 (derivs)
       real*8 derivs(3,*)
 c
 c
-c     distribute gradient on four-site water extra centers
+c     distribute the gradient on four-site water extra centers
 c
       if (use_rattle) then
          do i = 1, nwat4
diff --git a/source/readxyz.f b/source/readxyz.f
@@ -235,71 +235,73 @@ subroutine readxyz (ixyz)
 c
 c     for each atom, count and sort its attached atoms
 c
-      do i = 1, n
-         do j = maxval, 1, -1
-            if (i12(j,i) .ne. 0) then
-               n12(i) = j
-               goto 150
-            end if
+      if (.not. abort) then
+         do i = 1, n
+            do j = maxval, 1, -1
+               if (i12(j,i) .ne. 0) then
+                  n12(i) = j
+                  goto 150
+               end if
+            end do
+  150       continue
+            call sort (n12(i),i12(1,i))
          end do
-  150    continue
-         call sort (n12(i),i12(1,i))
-      end do
 c
 c     perform dynamic allocation of some local arrays
 c
-      nmax = 0
-      do i = 1, n
-         nmax = max(tag(i),nmax)
-         do j = 1, n12(i)
-            nmax = max(i12(j,i),nmax)
+         nmax = 0
+         do i = 1, n
+            nmax = max(tag(i),nmax)
+            do j = 1, n12(i)
+               nmax = max(i12(j,i),nmax)
+            end do
          end do
-      end do
-      allocate (list(nmax))
+         allocate (list(nmax))
 c
 c     check for scrambled atom order and attempt to renumber
 c
-      reorder = .false.
-      do i = 1, n
-         list(tag(i)) = i
-         if (tag(i) .ne. i)  reorder = .true.
-      end do
-      if (reorder) then
-         write (iout,160)
-  160    format (/,' READXYZ  --  Atom Labels not Sequential,',
-     &              ' Attempting to Renumber')
+         reorder = .false.
          do i = 1, n
-            tag(i) = i
-            do j = 1, n12(i)
-               i12(j,i) = list(i12(j,i))
-            end do
-            call sort (n12(i),i12(1,i))
+            list(tag(i)) = i
+            if (tag(i) .ne. i)  reorder = .true.
          end do
-      end if
+         if (reorder) then
+            write (iout,160)
+  160       format (/,' READXYZ  --  Atom Labels not Sequential,',
+     &                 ' Attempting to Renumber')
+            do i = 1, n
+               tag(i) = i
+               do j = 1, n12(i)
+                  i12(j,i) = list(i12(j,i))
+               end do
+               call sort (n12(i),i12(1,i))
+            end do
+         end if
 c
 c     perform deallocation of some local arrays
 c
-      deallocate (list)
+         deallocate (list)
 c
 c     check for atom pairs with identical coordinates
 c
-      clash = .false.
-      if (n .le. 10000)  call chkxyz (clash)
+         clash = .false.
+         if (n .le. 10000)  call chkxyz (clash)
 c
 c     make sure all atom connectivities are bidirectional
 c
-      do i = 1, n
-         do j = 1, n12(i)
-            k = i12(j,i)
-            do m = 1, n12(k)
-               if (i12(m,k) .eq. i)  goto 180
+         do i = 1, n
+            do j = 1, n12(i)
+               k = i12(j,i)
+               do m = 1, n12(k)
+                  if (i12(m,k) .eq. i)  goto 180
+               end do
+               write (iout,170)  k,i
+  170          format (/,' READXYZ  --  Check Connection of Atoms',
+     &                    i9,' and',i9)
+               call fatal
+  180          continue
             end do
-            write (iout,170)  k,i
-  170       format (/,' READXYZ  --  Check Connection of Atoms',
-     &                 i9,' and',i9)
-            call fatal
-  180       continue
          end do
-      end do
+      end if
       return
       end

Original file line number	Diff line number	Diff line change
`@@ -45,7 +45,7 @@ subroutine control`
`45`	`45`	`c check for control parameters on the command line`
`46`	`46`	`c`
`47`	`47`	`exist = .false.`
`48`		`- do i = 1, narg-1`
	`48`	`+ do i = 1, narg`
`49`	`49`	`string = arg(i)`
`50`	`50`	`call upcase (string)`
`51`	`51`	`if (string(1:2) .eq. '-D') then`
Original file line number	Diff line number	Diff line change
`@@ -150,7 +150,7 @@ subroutine rattle (dt,xold,yold,zold)`
`150`	`150`	`& ' Iterations')`
`151`	`151`	`end if`
`152`	`152`	`c`
`153`		`-c set four-site water extra centers from H-O-H coordinates`
	`153`	`+c set position and velocity of four-site water extra centers`
`154`	`154`	`c`
`155`	`155`	`do i = 1, nwat4`
`156`	`156`	`ia = iwat4(1,i)`
`@@ -162,6 +162,9 @@ subroutine rattle (dt,xold,yold,zold)`
`162`	`162`	`x(ia) = otermx(ib) + hterm(x(ic)+x(id))`
`163`	`163`	`y(ia) = otermy(ib) + hterm(y(ic)+y(id))`
`164`	`164`	`z(ia) = otermz(ib) + hterm(z(ic)+z(id))`
	`165`	`+ v(1,ia) = 0.0d0`
	`166`	`+ v(2,ia) = 0.0d0`
	`167`	`+ v(3,ia) = 0.0d0`
`165`	`168`	`end do`
`166`	`169`	`c`
`167`	`170`	`c apply group position and velocity constraints via exact reset`
`@@ -522,7 +525,7 @@ subroutine shake (xold,yold,zold)`
`522`	`525`	`& ' Iterations')`
`523`	`526`	`end if`
`524`	`527`	`c`
`525`		`-c set four-site water extra centers from H-O-H coordinates`
	`528`	`+c set the position of any four-site water extra centers`
`526`	`529`	`c`
`527`	`530`	`do i = 1, nwat4`
`528`	`531`	`ia = iwat4(1,i)`
`@@ -719,7 +722,7 @@ subroutine water4 (derivs)`
`719`	`722`	`real8 derivs(3,)`
`720`	`723`	`c`
`721`	`724`	`c`
`722`		`-c distribute gradient on four-site water extra centers`
	`725`	`+c distribute the gradient on four-site water extra centers`
`723`	`726`	`c`
`724`	`727`	`if (use_rattle) then`
`725`	`728`	`do i = 1, nwat4`