[gmx-users] Re: amber force field in Gromacs

[servaas]

> Dear Servaas,
> 
> In tleap did you really did:
> 
> TLEAP
> tleap -f leaprc.ff99SB
> ad = sequence { DA5 DA DA3 }
> saveamberparm da da_amber.top da_amber.crd
> 
> 
> If so, it's wrong, it should be:
> 
> saveamberparm ad da_amber.top da_amber.crd
>                         ^^^
> and not 'da'
Yes of course just a typo here
> 
> Besides, I tried to reproduce what you did using what I think would be
> fine and... everything went fine! Energies after minimisation in
> single and double were almost identical and trajectories diverted
> normally.
Did you also try running it in vacuum? In solvent the problems occur
only after a large number of steps. In vacuum they occur very fast. I
know vacuum is not the way to go, but I consider it as quick test, if a
short simulation in vacuum gives strange things it is an indication of a
problem with the parameters (force fiels or others...). Amber doesn't
give me any problems in vacuum, this gives me an indication that the
vacuum is not the problem here.
> 
> Please check what I did.
> 
> # begin commands
> 
> cat << EOF >| em.mdp
> define                   = -DFLEXIBLE
> integrator               = cg ; steep
> nsteps                   = 200
> constraints              = none
> emtol                    = 1000.0
> nstcgsteep               = 10 ; do a steep every 10 steps of cg
> emstep                   = 0.01 ; used with steep
> nstcomm                  = 1
> coulombtype              = PME
> ns_type                  = grid
> rlist                    = 1.0
> rcoulomb                 = 1.0
> rvdw                     = 1.4
> Tcoupl                   = no
> Pcoupl                   = no
> gen_vel                  = no
> nstxout                  = 0 ; write coords every # step
> optimize_fft             = yes
> EOF
> 
> 
> cat << EOF >| md.mdp
> integrator               = md
> nsteps                   = 1000
> dt                       = 0.002
> constraints              = all-bonds
> nstcomm                  = 1
> ns_type                  = grid
> rlist                    = 1.2
> rcoulomb                 = 1.1
> rvdw                     = 1.0
> vdwtype                  = shift
> rvdw-switch              = 0.9
> coulombtype              = PME-Switch
> Tcoupl                   = v-rescale
> tau_t                    = 0.1 0.1
> tc-grps                  = protein non-protein
> ref_t                    = 300 300
> Pcoupl                   = parrinello-rahman
> Pcoupltype               = isotropic
> tau_p                    = 0.5
> compressibility          = 4.5e-5
> ref_p                    = 1.0
> gen_vel                  = yes
> nstxout                  = 2 ; write coords every # step
> lincs-iter               = 2
> DispCorr                 = EnerPres
> optimize_fft             = yes
> EOF
> 
> 
> cat << EOF >| leap.in
> verbosity 1
> source leaprc.ff99SB
> ad = sequence { DA5 DA DA3 }
> solvatebox ad TIP3PBOX 10.0
> addions ad Na+ 5
> addions ad Cl- 3
> saveamberparm ad da_amber.top da_amber.crd
> savepdb ad DA.pdb
> quit
> EOF
> tleap -f leap.in >| leap.out
> 
> acpypi -x da_amber.crd -p da_amber.top -d
> 
> #Single precision
> grompp -f em.mdp -c da_amber_GMX.gro -p da_amber_GMX.top -o em.tpr
> mdrun -v -deffnm em
> 
> #Polak-Ribiere Conjugate Gradients converged to Fmax < 1000 in 22 steps
> #Potential Energy  = -6.2280516e+04
> #Maximum force     =  7.5868494e+02 on atom 98
> #Norm of force     =  1.0447179e+02
> 
> grompp -f md.mdp -c em.gro -p da_amber_GMX.top -o md.tpr
> mdrun -v -deffnm md
> 
> #Double precision
> grompp_d -f em.mdp -c da_amber_GMX.gro -p da_amber_GMX.top -o em.tpr
> mdrun_d -v -deffnm em
> 
> #Polak-Ribiere Conjugate Gradients converged to Fmax < 1000 in 22 steps
> #Potential Energy  = -6.22813514022256e+04
> #Maximum force     =  7.58238100790309e+02 on atom 98
> #Norm of force     =  1.04358667410458e+02
> 
> grompp_d -f md.mdp -c em.gro -p da_amber_GMX.top -o md.tpr
> mdrun_d -v -deffnm md
> 
> # end commands
> 
> Regards,
> 
> Alan