CHARMM Development Project Forums User Discussion & Questions Molecular Dynamics problems with long simulation time with drudes/VV2

Greetings all!
I am trying to calculate some dielectric properties of alcohols and am using drude oscillators with the VV2 integrator. I need to run the simulation for multiple nanoseconds. I can successfully run the MD for 350 ps (long enough to equilibrate), but the program aborts after times longer than that (which is what I need). If I restart after (or even shortly before) the MD aborts, the same error occurs.
The error is *ERROR IN SHAKE/ROLL* (see output below), and if I decrease the bomlev the error continues as:
DEBUG -- SHAKE/ROLL did NOT converge!
VV2 -- RMSDshake 1.643296882254187E-007
VV2 -- (Try increasing BTAU)
If I increase btau, then the volume grows to incredible sizes.
These errors occur while the Image and nonbond lists are updating. Could this be the source of error? Does anyone have experience with longer simulations that could help?
Any guidance is greatly appreciated!

Input:
++++++++++++++++++++++
READ SEQUence PRO1 128

GENErate SOLV first none last none SETUp warn DRUDE DMASS 0.4

OPEN READ UNIT 10 CARD NAME @coorfile1
READ COOR PDB UNIT 10

COOR SDRUDE
COOR COPY COMP

SHAKE BONH PARAM TOLEr 1.0E-10 NOFAST -
select ( resn PRO1 .and. .not. type d* ) end
COOR SHAKE

CRYStal DEFIne RECTangular 28 28 28 90. 90. 90.

CRYStal BUILd NOPER 0 CUTOff 16

NBONd -
inbfrq -1 imgfrq -1 ihbfrq 0 -
EWALD kappa 0.34 PMEWald fftx 24 ffty 24 fftz 24 order 6 spline -
ELEC atom cdiel fswitch -
VDW vatom vfswitch -
cutim 16 cutnb 16 ctofnb 14 ctonnb 12 EPS 1.0 -
lrc

MINI SD step 0.001 nstep 100 nprint 10

OPEN WRITE UNIT 29 CARD NAME ~/@dir/@idfile_amin.pdb
WRITE COOR PDB UNIT 29
CLOSE UNIT 29

tpcontrol nther 2 nstep 50 -
ther 1 tau 0.1 tref 289.15 select .not. type d* end -
ther 2 tau 0.005 tref 1 select type d* end -
baro btau 0.2 pref 1.00
!tpcontrol off

open unit 35 write form name ~/@dir/@idfile_@i.rst
open unit 32 write unfo name ~/@dir/@idfile_@i.trj
open unit 33 write form name ~/@dir/@idfile_@i.kunit
open unit 34 write unfo name ~/@dir/@idfile_@i.vel

DYNA vv2 start nstep @nstep1 time 0.001 ntrfrq 100 -
iprfrq 100 nprint 100 -
iasvel 1 firstt 298.15 finalt 298.15 -
iseed 314159 -
iunread -1 iunwrite 35 iuncrd 32 iunvel 34 kunit 33 -
nsavc 100 nsavv 100 isvfrq 100
stop
++++++++++++++++++++

Output error after 300 ps:
++++++++++++++++++++++
* * * AVERAGES FOR THE LAST 100 STEPS
AVER DYN: Step Time TOTEner TOTKe ENERgy TEMPerature
AVER PROP: GRMS HFCTote HFCKe EHFCor VIRKe
AVER INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
AVER EXTERN: VDWaals ELEC HBONds ASP USER
AVER IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElec
AVER EWALD: EWKSum EWSElf EWEXcl EWQCor EWUTil
AVER PRESS: VIRE VIRI PRESSE PRESSI VOLUme
---------- --------- --------- --------- --------- ---------
AVER> 100 379.35000 1864.15099 1072.01321 347.24302 300.35892
AVER PROP> 10.96833 1419.25623 0.00000 444.89476 1155.52612
AVER INTERN> 499.24864 516.51159 107.45168 458.24217 0.00000
AVER EXTERN> -141.36626 -263.47674 0.00000 0.00000 0.00000
AVER IMAGES> -347.77881 -721.73670 0.00000 0.00000 0.00000
AVER EWALD> 41.14430-218404.15549 218603.15865 0.00000 0.00000
AVER PRESS> -49.39759 -720.95316 222.63173 -27.66194 15277.93628
---------- --------- --------- --------- --------- ---------
Lattice Parameters> Averages for the last 100 steps:
Crystal Parameters : Crystal Type = RECT
AVER A = 24.81351 B = 24.81351 C = 24.81351
AVER Alpha = 90.00000 Beta = 90.00000 Gamma = 90.00000
AVER PIXX = 195.09 PIYY = -152.33 PIZZ = -125.75
AVER PIXY = -219.61 PIXZ = 27.22 PIYZ = 16.21
AVER Gradient Norm = 758.66566


* * * RMS FLUCTUATIONS FOR 100 STEPS
FLUC> 100 379.35000 0.22338 18.31053 15.63794 5.13019
FLUC PROP> 0.15080 8.26908 0.00000 8.21520 460.64601
FLUC INTERN> 20.61182 12.84091 3.75550 4.51328 0.00000
FLUC EXTERN> 10.40633 19.60115 0.00000 0.00000 0.00000
FLUC IMAGES> 4.85209 18.21459 0.00000 0.00000 0.00000
FLUC EWALD> 0.40814 0.00000 2.32724 0.00000 0.00000
FLUC PRESS> 298.77200 67.44139 1342.78597 285.64563 21.10045
---------- --------- --------- --------- --------- ---------
Lattice Parameters> RMS fluctuations for the last 100 steps:
Crystal Parameters : Crystal Type = RECT
FLUC A = 0.01142 B = 0.01142 C = 0.01142
FLUC Alpha = 0.00000 Beta = 0.00000 Gamma = 0.00000
FLUC PIXX = 946.96 PIYY = 1333.59 PIZZ = 1659.14
FLUC PIXY = 878.57 PIXZ = 632.40 PIYZ = 817.23
FLUC Gradient Norm = 690.21134


DRIFT/STEP (LAST-TOTAL): -1.57326775E-03 -1.57326775E-03
EPROP(EPOT) AT STEP 0 : 1864.2304 1864.2304
CORR. COEFFICIENT : -0.20330724 -0.20330724
UPDECI: Image update at step 379106
UPDECI: Nonbond update at step 379106
UPDECI: Image update at step 379124
UPDECI: Nonbond update at step 379124
UPDECI: Image update at step 379125
UPDECI: Nonbond update at step 379125
UPDECI: Image update at step 379126
UPDECI: Nonbond update at step 379126
UPDECI: Image update at step 379127
UPDECI: Nonbond update at step 379127
** ERROR IN ShakeRoll ** DEVIATION IN SHAKE TOO LARGE
NITER= 3 LL= 267 I= 597 J= 599
TOLER= 0.9409000000 DIFF= 7.2289133683
RRPR= -1.8533666500
X(I) -5.5339163901 4.3065370393 -10.0417178202
XREF(I) -5.5300331959 4.3598011072 -10.0748854954

Since the edge lengths are all the same, I suggest using CUBIC instead of RECTangular for the CRYStal lattice type.

Technically, it's an error to specify both EWALD and FSWITCH, but CHARMM assumes you meant EWALD in this case.

However, I've had similar ShakeRoll problems with VV2/DRUDE, which I was unable to resolve, for some PEG/alkane lipids.

I thought RECTangular was the CUBIC equivalent for VV2. (compatibility issue) I will try it and see.

Thanks!

No, RECTangular was added as a special ORTHorhombic case with reduced degrees of freedom, such that ratio of edge lengths is maintained. With ORTHorhombic, one edge can shrink as another one grows, allowing shape changes.

I'm quite surprised that increasing BTAU would cause your system to blow up. BTAU is the characteristic response time of the barostat, and a longer BTAU means a slower "piston". (Actually, if you make BTAU large enough, say, 1.0E8, the volume of your box shouldn't change at all. Quite the opposite of what you see...)

Do you get the same results/problems with a CUBIC crystal?

Sorry for the delay, the flu got me (regular, not the swine)

I replaced RECTangular with CUBIc:
RECT simulation aborts at 209 ps; CUBIC simulation aborts at 46 ps. No other variables were changed. There is also no indication as to why the simulation terminated. I do, however, get a better volume when using CUBIC (proper density).

When I increase BTAU, I adjust it from 0.1 to 0.5. At 0.5 it blows up, but I have not tried higher than that. Is it appropriate to use such a large BTAU, or does that negate the point of even having a separate barostat?

Thanks for all of your guidance. (My advisor thanks you as well)