[gmx-users] Free energy perturbation parameters sc-alpha sc-power sc-sigma
Francesca Vitalini
francesca.vitalini11 at gmail.com
Thu Jan 23 10:14:37 CET 2014
Hi,
thank you for the response,
Il giorno 23/gen/2014, alle ore 06:35, bipin singh <bipinelmat at gmail.com> ha scritto:
>> I am trying to perform a series of simulations where I gradually transform
>> a valine capped amino acid (Ac-VAL-NHMe) into an alanine capped amino acid
>> (Ac-ALA-NHMe), ie I have to turn on the HG and transform CG into HB. I have
>> defined 6 lambda values at which I would like to equilibrate my system
>> (lambda=0,0.2,0.40.6,0.8,1)
>
> Would it not be better if you turn off the vdw and coul. terms of CG (1
> atom) and HG (two atoms) and transform one of the HG to HB while doing the
> transformation from valine to alanine, rather than mutating the CG to HB.
>
If I transform a HG into a HB then won't I have problems with the bonds definitions? How should I change my topology in order to define that in the end state the mutated HG has to be bound to CB instead?
Thank you very much.
Francesca
>
> On Wed, Jan 22, 2014 at 9:40 PM, Francesca Vitalini <
> francesca.vitalini11 at gmail.com> wrote:
>
>> Dear Gromacs users,
>>
>> I am trying to perform a series of simulations where I gradually transform
>> a valine capped amino acid (Ac-VAL-NHMe) into an alanine capped amino acid
>> (Ac-ALA-NHMe), ie I have to turn on the HG and transform CG into HB. I have
>> defined 6 lambda values at which I would like to equilibrate my system
>> (lambda=0,0.2,0.40.6,0.8,1)
>>
>> I am using GROMACS.4.5.5. and OPLS-AA force field.
>>
>> I have been following the tutorial by Justin Lemkul
>>
>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html
>>
>> I have problems with the nvt equilibration with lambda=0.8 and lambda=1 as
>> my systems systematically crashes. I believe that the reason is in the
>> tuning of the parameters sc-alpha sc-power sc-sigma in the mdp file.
>> I understand from the manual that sc_sigma should not be the default option
>> as I have 6 hydrogens disappearing, and consequently the sc-alpha sc-power
>> parameters will be effected as the interactions of the hydrogens in the
>> soft core state will get stronger.
>> However I have no experience in how to select these parameters and would
>> like to ask, if anyone has more experience, which could be plausible
>> values.
>>
>> Thank you very much,
>>
>> Francesca
>>
>>
>> For completeness I attach here a copy of my topology and mdp files.
>>
>>
>> TOPOLOGY
>>
>>
>> ; Include forcefield parameters
>> #include "oplsaa.ff/forcefield.itp"
>>
>> [ moleculetype ]
>> ; Name nrexcl
>> Protein_chain_A 3
>>
>> [ atoms ]
>> ; nr type resnr residue atom cgnr charge mass
>> typeB chargeB massB
>> ; residue 1 ACE rtp ACE q 0.0
>> 1 opls_135 1 ACE CH3 1 -0.18 12.011
>> opls_135 -0.18 12.011
>> 2 opls_140 1 ACE HH31 1 0.06 1.008
>> opls_140 0.06 1.008
>> 3 opls_140 1 ACE HH32 1 0.06 1.008
>> opls_140 0.06 1.008
>> 4 opls_140 1 ACE HH33 1 0.06 1.008
>> opls_140 0.06 1.008
>> 5 opls_235 1 ACE C 2 0.5 12.011
>> opls_235 0.5 12.011
>> 6 opls_236 1 ACE O 2 -0.5 15.9994
>> opls_236 -0.5 15.9994
>> ; residue 2 VAL rtp VAL q 0.0
>> 7 opls_238 2 VAL N 3 -0.5 14.0067
>> opls_238 -0.5 14.0067
>> 8 opls_241 2 VAL H 3 0.3 1.008
>> opls_241 0.3 1.008
>> 9 opls_224B 2 VAL CA 3 0.14 12.011
>> opls_224B 0.14 12.011
>> 10 opls_140 2 VAL HA 3 0.06 1.008
>> opls_140 0.06 1.008
>> 11 opls_137 2 VAL CB 4 -0.06 12.011
>> opls_137 -0.18 12.011
>> 12 opls_140 2 VAL HB 4 0.06 1.008
>> opls_140 0.06 1.008
>> 13 opls_135 2 VAL CG1 5 -0.18 12.011
>> opls_140 0.06 1.008
>> 14 opls_140 2 VAL HG11 5 0.06 1.008
>> opls_140 0.00 0.000
>> 15 opls_140 2 VAL HG12 5 0.06 1.008
>> opls_140 0.00 0.000
>> 16 opls_140 2 VAL HG13 5 0.06 1.008
>> opls_140 0.00 0.000
>> 17 opls_135 2 VAL CG2 6 -0.18 12.011
>> opls_140 0.06 1.008
>> 18 opls_140 2 VAL HG21 6 0.06 1.008
>> opls_140 0.00 0.000
>> 19 opls_140 2 VAL HG22 6 0.06 1.008
>> opls_140 0.00 0.000
>> 20 opls_140 2 VAL HG23 6 0.06 1.008
>> opls_140 0.00 0.000
>> 21 opls_235 2 VAL C 7 0.5 12.011
>> opls_235 0.5 12.011
>> 22 opls_236 2 VAL O 7 -0.5 15.9994
>> opls_236 -0.5 15.9994
>> ; residue 3 NAC rtp NAC q 0.0
>> 23 opls_238 3 NAC N 8 -0.5 14.0067
>> opls_238 -0.5 14.0067
>> 24 opls_241 3 NAC H 8 0.3 1.008
>> opls_241 0.3 1.008
>> 25 opls_242 3 NAC CH3 9 0.02 12.011
>> opls_242 0.02 12.011
>> 26 opls_140 3 NAC HH31 9 0.06 1.008
>> opls_140 0.06 1.008
>> 27 opls_140 3 NAC HH32 9 0.06 1.008
>> opls_140 0.06 1.008
>> 28 opls_140 3 NAC HH33 9 0.06 1.008
>> opls_140 0.06 1.008
>>
>> [ bonds ]
>> ; ai aj funct c0 c1 c2 c3
>> 1 2 1
>> 1 3 1
>> 1 4 1
>> 1 5 1
>> 5 6 1
>> 5 7 1
>> 7 8 1
>> 7 9 1
>> 9 10 1
>> 9 11 1
>> 9 21 1
>> 11 12 1
>> 11 13 1
>> 11 17 1
>> 13 14 1
>> 13 15 1
>> 13 16 1
>> 17 18 1
>> 17 19 1
>> 17 20 1
>> 21 22 1
>> 21 23 1
>> 23 24 1
>> 23 25 1
>> 25 26 1
>> 25 27 1
>> 25 28 1
>>
>> [ pairs ]
>> ; ai aj funct c0 c1 c2 c3
>> 1 8 1
>> 1 9 1
>> 2 6 1
>> 2 7 1
>> 3 6 1
>> 3 7 1
>> 4 6 1
>> 4 7 1
>> 5 10 1
>> 5 11 1
>> 5 21 1
>> 6 8 1
>> 6 9 1
>> 7 12 1
>> 7 13 1
>> 7 17 1
>> 7 22 1
>> 7 23 1
>> 8 10 1
>> 8 11 1
>> 8 21 1
>> 9 14 1
>> 9 15 1
>> 9 16 1
>> 9 18 1
>> 9 19 1
>> 9 20 1
>> 9 24 1
>> 9 25 1
>> 10 12 1
>> 10 13 1
>> 10 17 1
>> 10 22 1
>> 10 23 1
>> 11 22 1
>> 11 23 1
>> 12 14 1
>> 12 15 1
>> 12 16 1
>> 12 18 1
>> 12 19 1
>> 12 20 1
>> 12 21 1
>> 13 18 1
>> 13 19 1
>> 13 20 1
>> 13 21 1
>> 14 17 1
>> 15 17 1
>> 16 17 1
>> 17 21 1
>> 21 26 1
>> 21 27 1
>> 21 28 1
>> 22 24 1
>> 22 25 1
>> 24 26 1
>> 24 27 1
>> 24 28 1
>>
>> [ angles ]
>> ; ai aj ak funct c0 c1
>> c2 c3
>> 2 1 3 1
>> 2 1 4 1
>> 2 1 5 1
>> 3 1 4 1
>> 3 1 5 1
>> 4 1 5 1
>> 1 5 6 1
>> 1 5 7 1
>> 6 5 7 1
>> 5 7 8 1
>> 5 7 9 1
>> 8 7 9 1
>> 7 9 10 1
>> 7 9 11 1
>> 7 9 21 1
>> 10 9 11 1
>> 10 9 21 1
>> 11 9 21 1
>> 9 11 12 1
>> 9 11 13 1
>> 9 11 17 1
>> 12 11 13 1
>> 12 11 17 1
>> 13 11 17 1
>> 11 13 14 1
>> 11 13 15 1
>> 11 13 16 1
>> 14 13 15 1
>> 14 13 16 1
>> 15 13 16 1
>> 11 17 18 1
>> 11 17 19 1
>> 11 17 20 1
>> 18 17 19 1
>> 18 17 20 1
>> 19 17 20 1
>> 9 21 22 1
>> 9 21 23 1
>> 22 21 23 1
>> 21 23 24 1
>> 21 23 25 1
>> 24 23 25 1
>> 23 25 26 1
>> 23 25 27 1
>> 23 25 28 1
>> 26 25 27 1
>> 26 25 28 1
>> 27 25 28 1
>>
>> [ dihedrals ]
>> ; ai aj ak al funct c0 c1
>> c2 c3 c4 c5
>> 2 1 5 6 3
>> 2 1 5 7 3
>> 3 1 5 6 3
>> 3 1 5 7 3
>> 4 1 5 6 3
>> 4 1 5 7 3
>> 1 5 7 8 3
>> 1 5 7 9 3
>> 6 5 7 8 3
>> 6 5 7 9 3
>> 5 7 9 10 3
>> 5 7 9 11 3
>> 5 7 9 21 3
>> 8 7 9 10 3
>> 8 7 9 11 3
>> 8 7 9 21 3
>> 7 9 11 13 3 dih_VAL_chi1_N_C_C_C
>> 7 9 11 17 3 dih_VAL_chi1_N_C_C_C
>> 21 9 11 13 3 dih_VAL_chi1_C_C_C_CO
>> 21 9 11 17 3 dih_VAL_chi1_C_C_C_CO
>> 7 9 11 12 3
>> 10 9 11 12 3
>> 10 9 11 13 3
>> 10 9 11 17 3
>> 21 9 11 12 3
>> 7 9 21 22 3
>> 7 9 21 23 3
>> 10 9 21 22 3
>> 10 9 21 23 3
>> 11 9 21 22 3
>> 11 9 21 23 3
>> 9 11 13 14 3
>> 9 11 13 15 3
>> 9 11 13 16 3
>> 12 11 13 14 3
>> 12 11 13 15 3
>> 12 11 13 16 3
>> 17 11 13 14 3
>> 17 11 13 15 3
>> 17 11 13 16 3
>> 9 11 17 18 3
>> 9 11 17 19 3
>> 9 11 17 20 3
>> 12 11 17 18 3
>> 12 11 17 19 3
>> 12 11 17 20 3
>> 13 11 17 18 3
>> 13 11 17 19 3
>> 13 11 17 20 3
>> 9 21 23 24 3
>> 9 21 23 25 3
>> 22 21 23 24 3
>> 22 21 23 25 3
>> 21 23 25 26 3
>> 21 23 25 27 3
>> 21 23 25 28 3
>> 24 23 25 26 3
>> 24 23 25 27 3
>> 24 23 25 28 3
>>
>> [ dihedrals ]
>> ; ai aj ak al funct c0 c1
>> c2 c3
>> 1 7 5 6 1 improper_O_C_X_Y
>> 5 9 7 8 1 improper_Z_N_X_Y
>> 9 23 21 22 1 improper_O_C_X_Y
>> 21 25 23 24 1 improper_Z_N_X_Y
>>
>> ; Include Position restraint file
>> #ifdef POSRES
>> #include "Ac_V_NHMe/Ac_V_NHMe.itp"
>> #endif
>>
>> ; Include water topology
>> #include "oplsaa.ff/tip3p.itp"
>>
>> #ifdef POSRES_WATER
>> ; Position restraint for each water oxygen
>> [ position_restraints ]
>> ; i funct fcx fcy fcz
>> 1 1 1000 1000 1000
>> #endif
>>
>> ; Include topology for ions
>> #include "oplsaa.ff/ions.itp"
>>
>> [ system ]
>> ; Name
>> ACETYL-VALINE-METHYLAMIDE in water
>>
>> [ molecules ]
>> ; Compound #mols
>> Protein_chain_A 1
>> SOL 682
>>
>>
>>
>>
>> NVT_VDW-Q.mdp
>> here I intend to switch off the coulomb interactions.
>>
>> ; Run control
>> integrator = sd ; Langevin dynamics
>> tinit = 0
>> dt = 0.002
>> nsteps = 10000 ; 20 ps
>> nstcomm = 100
>> ; Output control
>> nstxout = 10000
>> nstvout = 10000
>> nstfout = 0
>> nstlog = 1000
>> nstenergy = 500
>> nstxtcout = 500
>> xtc-precision = 500
>> xtc_grps = Protein
>> ; Neighborsearching and short-range nonbonded interactions
>> nstlist = 5
>> ns_type = grid
>> pbc = xyz
>> rlist = 1.0
>> ; Electrostatics
>> coulombtype = PME
>> rcoulomb = 1.0
>> ; van der Waals
>> vdw-type = switch
>> rvdw-switch = 0.8
>> rvdw = 0.9
>> ; Apply long range dispersion corrections for Energy and Pressure
>> DispCorr = EnerPres
>> ; Spacing for the PME/PPPM FFT grid
>> fourierspacing = 0.12
>> ; EWALD/PME/PPPM parameters
>> pme_order = 6
>> ewald_rtol = 1e-06
>> epsilon_surface = 0
>> optimize_fft = no
>> ; Temperature coupling
>> ; tcoupl is implicitly handled by the sd integrator
>> tc_grps = system
>> tau_t = 1.0
>> ref_t = 300
>> ; Pressure coupling is off for NVT
>> Pcoupl = No
>> tau_p = 0.5
>> compressibility = 4.5e-05
>> ref_p = 1.0
>> ; Free energy control stuff
>> free_energy = yes
>> init_lambda = 0.8
>> delta_lambda = 0
>> sc-alpha = 0.5
>> sc-power = 1.0
>> sc-sigma = 0.005
>> couple-moltype = Protein_chain_A ; name of moleculetype to
>> decouple
>> couple-lambda0 = vdw-q ; only van der Waals interactions
>> couple-lambda1 = vdw ; turn off everything, in this case only
>> vdW
>> couple-intramol = no
>> ; Generate velocities to start
>> gen_vel = yes
>> gen_temp = 300
>> gen_seed = -1
>> ; options for bonds
>> constraints = h-bonds ; we only have C-H bonds here
>> ; Type of constraint algorithm
>> constraint-algorithm = lincs
>> ; Do not constrain the starting configuration
>> continuation = no
>> ; Highest order in the expansion of the constraint coupling matrix
>> lincs-order = 12
>>
>>
>> NVT_VDW.mdp
>> here I intend to switch off the Van der Waals interactions. The only
>> difference with the previous file is in the Free energy parameters, so I
>> will report only those.
>>
>> ; Free energy control stuff
>> free_energy = yes
>> init_lambda = 0.8
>> delta_lambda = 0
>> sc-alpha = 0.5
>> sc-power = 1.0
>> sc-sigma = 0.3
>> couple-moltype = Protein_chain_A ; name of moleculetype to
>> decouple
>> couple-lambda0 = vdw ; only van der Waals interactions
>> couple-lambda1 = none ; turn off everything, in this case
>> only vdW
>> --
>> 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.
>>
>
>
>
> --
>
>
>
> *--------------------Thanks and Regards,Bipin Singh*
> --
> 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