[gmx-users] problems with reaction field

[hugo verli]

Hi to all,

I will return to previously discussed subjects on this list. I have simulated
a system composed by protein, ligand, ions, and spc water with PME and
berendsen pressure scaling for 5ns and everything works fine. To facilitate
the use of LIE approximation for free energy calculation (I intent to compare
the values obtained for several molecules) I decided to run another simulation
of the same system using reaction field and pressure scaling. Then some
problems appeared. Initially the error

> Warning: pressure scaling more than 1%, mu: 1.04684 1.04684=20

occurred. Searching the GROMACS mail list I found a recomendation  to test
different values of tau_p. So I modificate this value from 1.0, 5.0, 10.0 and
finally to 20.0, and the error still occur. Then I saw another recomendation
in the mail list, to run first a NVT simulation and if everything is fine run
a NPT simulation. A NVT simulation of about 100ps works fine, but it crashes
with the same error when I introduce the pressure scaling. So I increase the
NVT simulation to about 500ps. The simulation do not crashe but display the error

> Constraint error in algorithm Lincs at step 159754
> Wrote pdb files with previous and current coordinates
> Large VCM(group rest):  10054.12793, -14113.24219, -12220.23730, ekin-cm: 
> 81059018440704.00000

And the system appeared to have exploded. Also some parameters at the .log
file, as well as the coordinates  in .gro file are substituted by a "nan"
character:

>   Energies (kJ/mol)
>          Angle    Proper Dih.  Improper Dih.          LJ-14     Coulomb-14
>            nan            nan            nan            nan            nan
>        LJ (SR)   Coulomb (SR)   Coulomb (LR)      Potential    Kinetic En.
>    0.00000e+00    0.00000e+00    0.00000e+00            nan            nan
>   Total Energy    Temperature Pressure (bar)
>            nan            nan            nan
>
>           Step           Time         Lambda      Annealing
>         234800      499.60004        0.00000        1.00000
>
>   Rel. Constraint Deviation:  Max    between atoms     RMS
>       Before LINCS         0.000000      1      2   nan
>        After LINCS         0.000000      1      2   nan

The .mdp file used was:

> title               =  Yo
> cpp                 =  /lib/cpp
> define              =  -DFLEX_SPC
> constraints         =  all-bonds
> integrator          =  md
> tinit               =  30.0
> dt                  =  0.002    ; ps !
> nsteps              =  235000   ; total 30-500 ps.
> nstcomm             =  1
> nstxout             =  250
> nstvout             =  1000
> nstfout             =  0
> nstlog              =  100
> nstenergy           =  100
> nstlist             =  10
> ns_type             =  grid
> coulombtype         =  Generalized-Reaction-Field
> rlist               =  0.9
> rcoulomb            =  1.5
> rvdw                =  0.9
> fourierspacing      =  0.12
> optimize_fft        =  no
> epsilon_r           =  80
> ; Berendsen temperature coupling is on in four groups
> Tcoupl              =  berendsen
> tc-grps             =  SOL Cl Protein DRG
> tau_t               =  0.1 0.1 0.1 0.1
> ref_t               =  310 310 310 310
> ; Energy monitoring
> energygrps          =  SOL Cl Protein DRG
> ; Isotropic pressure coupling is now on
> Pcoupl              =  no
> Pcoupltype          = isotropic
> tau_p               =  5.0
> compressibility     =  4.5e-5
> ref_p               =  1.0
> ; Generate velocites is off at 100 K.
> gen_vel             =  no
> gen_temp            =  310.0
> gen_seed            =  173529

Did I made something wrong? The topology obtained for the ligand is ok, and
the simulation runs without any problem in PME. Curiously I have obtained the
same errors with distinct protein-protein system. 

I was thinking to run a reaction field simulation using the last structure
obtained from PME simulation. Does this make sense? 

Many thanks in advance form any help,

MSc. Hugo Verli
Centro de Biotecnologia
Universidade Federal do Rio Grande do Sul