********> bugfix.77 Author: Tom Cheatham and Bill Ross Date: 9/24/96 Programs: Sander Severity: (As much enhancement as bugfix) Problem: Energy drains during dynamics, coupled with Berendsen temperature coupling can lead to a slow growth in the center of mass velocity, resulting in a 'flying block of ice' as this form of kinetic energy predominates (this has been observed in Ewald simulations). This fix enables center-of-mass motion removal when using periodic boundaries, via the NTCM and NSCM input variables (these settings only affected vacuum systems before). Fix: Make the following change to sander.f: ------------------------------------------------------------------------- *** OLD sander.f --- NEW sander.f *************** *** 411,416 **** --- 411,450 ---- + X(L40),X(L35),X(L20),TEMPI,HEAT,DT,INIT,IG,iscale, + scalm) IF(BELLY) CALL BELLYF(NATOM,IX(I62),X(L40)) + c + c --- remove com motion if req'd --- + c + if (imin.eq.0 .and. ntcm.ne.0 .and. ntb.ne.0) then + vx = 0.d0 + vy = 0.d0 + vz = 0.d0 + + j = 1 + tmass = 0.0d0 + do i = 1, 3*natom,3 + amass = 1.0d0 / x(l20+j-1) + tmass = tmass + amass + + vx = vx + amass * x(l40+i-1) + vy = vy + amass * x(l40+i) + vz = vz + amass * x(l40+i+1) + j = j + 1 + enddo + + vx = vx / tmass + vy = vy / tmass + vz = vz / tmass + + vel = sqrt(vx*vx + vy*vy + vz*vz) + if ( vel.gt.0.01d0 ) + . write (6,'(a,f9.2)') ' Removed COM velocity: ', vel + + do i = 1, 3*natom, 3 + x(l40+i-1) = x(l40+i-1) - vx + x(l40+i) = x(l40+i) - vy + x(l40+i+1) = x(l40+i+1) - vz + enddo + endif C C If we're reading NMR restraints/weight changes, read them now: C ------------------------------------------------------------------------- Fix: Make the following change to runmd.f: ------------------------------------------------------------------------- *** OLD runmd.f --- NEW runmd.f *************** *** 1040,1045 **** --- 1040,1091 ---- IF(BELLY) CALL BELLYF(NR,IX(I62),V) c c------------------------------------------------------------------- + c Step 5a: zero COM velocity if requested; used for preventing + c ewald 'block of ice flying thru space' phenomenon + c------------------------------------------------------------------- + c + if (mod(nstep+1,nscm).eq.0 .and. NTB.ne.0) then + c + vcmx = 0.d0 + vcmy = 0.d0 + vcmz = 0.d0 + + j = 1 + do i = 1, 3*natom,3 + amass = 1.0d0 / winv(j) + + vcmx = vcmx + amass * v(i) + vcmy = vcmy + amass * v(i+1) + vcmz = vcmz + amass * v(i+2) + + j = j + 1 + enddo + + vcmx = vcmx / tmas + vcmy = vcmy / tmas + vcmz = vcmz / tmas + + vel2 = vcmx*vcmx + vcmy*vcmy + vcmz*vcmz + atempdrop = 0.5d0 * tmas * vel2 / fac(1) + vel = sqrt(vel2) + #ifdef MPI + if ( master .and. vel.gt.0.01d0 ) then + #else + if ( vel.gt.0.01d0 ) then + #endif + write (6,'(a,f9.2,a,f9.2,a)') + . ' Removed COM velocity ',vel,' approx temp drop ', + . atempdrop, ' K' + endif + + do i = 1, 3*natom, 3 + v(i) = v(i) - vcmx + v(i+1) = v(i+1) - vcmy + v(i+2) = v(i+2) - vcmz + enddo + endif + c + c------------------------------------------------------------------- c Step 6: scale coordinates if constant pressure run: c------------------------------------------------------------------- c ------------------------------------------------------------------------- Temporary workarounds: NSNB (pairlist update frequency) = 1, TOL (SHAKE tolerance) = 0.000001, constant volume. --