[gmx-users] maintaining stability when uncoupling ion pairs
Michael Shirts
mrshirts at gmail.com
Thu Jul 9 23:35:07 CEST 2015
I would probably do this not using the couple-moltype functionality,
since you want to change TWO molecules. GROMACS thinks they are part
of the same molecule, and thus they should not move too far away (I'm
guessing). Instead, you can do this through the definition of A and B
types (see documentation)
On Thu, Jul 9, 2015 at 2:19 PM, andrew biedermann <amb4ht at virginia.edu> wrote:
> Hello,
>
> I want to calculate the free energy difference associated with decoupling a
> NaCl pair from aqueous solution, but am running into a strange error. I
> believe the issue is in how I define the molecule type in the topology
> file, since the simulation is stable when removing individual ions, or
> during regular production runs with the same mdp settings. I get the
> following error soon after submitting the job (5000 steps in):
>
> A list of missing interactions:
> LJC Pairs NB of 1 missing 1
> exclusions of 16669 missing 1
>
> Molecule type 'tNC'
> the first 10 missing interactions, except for exclusions:
> LJC Pairs NB atoms 1 2 global 11113 11114
>
> -------------------------------------------------------
> Program mdrun, VERSION 5.0.4
> Source code file:
> /nv/blue/amb4ht/downloads/gromacs-5.0.4/src/gromacs/mdlib/domdec_top.c,
> line: 394
>
> Fatal error:
> 2 of the 19448 bonded interactions could not be calculated because some
> atoms involved moved further apart than the multi-body cut-off distance
> (1.74968 nm) or the two-body cut-off distance (1.74968 nm), see option
> -rdd, for pairs and tabulated bonds also see option -ddcheck
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
> -------------------------------------------------------
>
> I have modified my topology file with the following lines:
>
> [ moleculetype ]
> ; molname nrexcl
> tNC 1
>
> [ atoms ]
> ; id at type res nr residu name at name cg nr charge
> 1 Na 1 NA NA 1 1.00000
> 2 Cl 1 CL CL 2 -1.00000
>
> [ bonds ]
>
> [ angles ]
>
> And my mdp file is:
>
> ; Run control
> integrator = sd ; Langevin dynamics
> tinit = 0
> dt = 0.002
> nsteps = 5000000 ; 10 ns
> ; Output control
> nstxout = 2500
> nstvout = 2500
> nstfout = 0
> nstlog = 2500
> nstenergy = 2500
> nstxout-compressed = 0
> Bond parameters
> constraint_algorithm = lincs ; holonomic constraints
> ;constraints = h-bonds ; all bonds (even heavy atom-H
> bonds) constrained
> ;MRS: use h-bonds
> lincs_iter = 2 ; accuracy of LINCS
> lincs_order = 12 ; also related to accuracy
> lincs_warnangle = 30
> continuation = yes
> ; Neighborsearching and short-range nonbonded interactions
> ; Neighborsearching
> cutoff-scheme = Group ; Group works with free energy
> rlist = 1.2 ; cut-off distance for the short-range neighbor list
> ns_type = grid ; search neighboring grid cells
> nstlist = 10
> ; Electrostatics
> coulombtype = PME
> rcoulomb = 1.2
> ; van der Waals
> vdw-type = Cutoff
> vdw-modifier = Potential-switch
> rvdw-switch = 0.90
> rvdw = 0.95
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr = EnerPres
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing = 0.10
> ; EWALD/PME/PPPM parameters
> pme_order = 6
> ewald_rtol = 1e-06
> epsilon_surface = 0
> ; Temperature coupling
> tc-grps = System
> tau_t = 2.0 ; time constant, in ps
> ref_t = 298 ; reference temperature in K
> ; Pressure coupling is on for NPT
> Pcoupl = Parrinello-Rahman
> tau_p = 2.0
> compressibility = 4.5e-05
> ref_p = 1.0
> ;############### Free energy control stuff ##############################
> free_energy = yes
> init_lambda_state = 0
> delta_lambda = 0
> calc_lambda_neighbors = 1 ; only immediate neighboring windows
> ; Vectors of lambda specified here
> ; Each combination is an index that is retrieved from init_lambda_state for
> each simulation
> ; init_lambda_state 0 1 2 3 4 5 6 7 8
> 9 10 11 12 13 14 15 16 17 18 19 20
> vdw_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
> coul_lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80
> 0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
> ; We are not transforming any bonded or restrained interactions
> bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> ; Masses are not changing (particle identities are the same at lambda = 0
> and lambda = 1)
> mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> ; Not doing simulated temperting here
> temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> ; Options for the decoupling
> sc-alpha = 0.5
> sc-coul = no ; linear interpolation of Coulomb
> sc-power = 1.0
> sc-sigma = 0.3
> couple-moltype = tNC ; name of moleculetype to decouple
> couple-lambda0 = vdw-q ; all interactions at lambda 0
> couple-lambda1 = none ; turn off everything at lambda 1
> couple-intramol = no
> nstdhdl = 100
>
> Thanks for the help,
>
> Drew
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
More information about the gromacs.org_gmx-users
mailing list