CHARMM Development Project Forums User Discussion & Questions Molecular Dynamics Problems with MD using polarizable ff

First, Thanks to all of you in this forum. My last 6 months has been so much easier reading all or your problems/solutions. Thanks a million!
Now, my problem:
I have been attempting to use a polarizable force field (drude oscillators) to model a box of propanol with dissociated Li and Cl. After many failed attempts, I keep getting *Error in ShakeRoll* which I have read is more of a symptom to a problem than the actual problem. I am using the vv2 integrator with the temperature control and two different thermostats (one for the propanol/li/cl and one for the drude oscillators). Has anyone had any experience using the drude oscillators and had similar problems? Any guidance is appreciated. (I can no longer bang my head against the wall) Below is the important part of my input file starting after reading in the rtf,etc.
Many Thanks - Allison

GENErate PROP first none last none SETUp warn DRUDE DMASS 0.4

OPEN READ UNIT 10 CARD NAME ~/generate/propanol_licl_128.pdb
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 45 45 45 90. 90. 90.

CRYStal BUILd NOPER 0 CUTOff 14.0

ENERgy -
inbfrq -1 imgfrq -1 ihbfrq 0 -
ewald pmewald kappa 0.34 fftx 54 ffty 54 fftz 54 order 6 -
switch atom vatom vfswitch -
cutnb 13.0 ctofnb 12.0 ctonnb 8.0 qcor 0.0

! Turn on image centering of solvent
image byres xcen 0.0 ycen 0.0 zcen 0.0 -
sele resn PRO1 .and. resn cla .and. resn lit END

mini sd step 0.001 nstep 600 nprint 50

OPEN WRITE UNIT 29 CARD NAME "~/propanol/prop3_amin.pdb"
WRITE COOR PDB UNIT 29
CLOSE UNIT 29

tpcontrol nther 2 cmdamp 10.0 nstep 20 -
ther 1 tau 0.1 tref 298.15 select .not. type D* end -
ther 2 tau 0.005 tref 0.1 select type D* end -
baro btau 0.2 pref 1.00

open unit 31 write form name ~/propanol/prop3_0.rst
open unit 32 write unfo name ~/propanol/prop3_0.trj
open unit 33 write form name ~/propanol/prop3_0.ene

!MD simulation
DYNA vv2 start nstep 1000 time 0.001 ntrfrq 100 -
iprfrq -1 nprint 100 -
iasvel 1 firstt 298.15 finalt 298.15 -
iseed 123456 -
iunread -1 iunwrite 31 iuncrd 32 iunvel -1 kunit 33 -
nsavc 1000 nsavv -1 isvfrq 1000

coor stat select all end

STOP

Your input looks OK.

I'll presume your initial coordinates (~/generate/propanol_licl_128.pdb) are good. Are they still OK after the minimization?

Some random tricks to get the simulation going:
- Try a different minimizer (ABNR instead of SD)
- Use a higher temperature for the Drudes (say, 1K instead of 0.1K)
- Use a very short time step (say, 0.1fs instead of 1fs)

If you can do a few hundreds steps using one of those tricks, save the coordinates and start the MD over with your regular settings. The goal of this first MD run is just to get rid of some possible bad conformations in your system.

Please write back to the forum if you still can't get the system running -- or to tell us what has worked.

Thanks!
I used ABNR rather than SD for the minimization and used a timestep of .1 fs for the initial dynamics simulation. Then, with the restart, I was able to successfully have a timestep of 1fs.
Everything seems to be working for now. Thanks again for your help.
Allison