[gmx-users] Free energy perturbation parameters sc-alpha sc-power sc-sigma
bipin singh
bipinelmat at gmail.com
Thu Jan 23 10:56:41 CET 2014
Ok, now I understand the transformation scheme. I think for doing the
mentioned transformation you have to do the following things:
1) In the first set of transformation, transform the HG1 (1,2 and 3) and
HG2 (1,2 and 3) to dummy atoms (you have written "I have to turn on the HG"
can you explain why?) by first relaxing the coulombic terms (first
transformation) followed by another transformation to relax the vdw terms
(in a separate transformation).
2) Then transform your CG1 and CG2 atoms to HB atoms in another set of
transformations.
During these transformation you have to choose the spacing between the
lambda values wisely based on some previous studies reported in literature
for similar type of transformation. Also the soft core parameters during
the (1) and (2) transformation should be set accrodingly.
On Thu, Jan 23, 2014 at 2:44 PM, Francesca Vitalini <
francesca.vitalini11 at gmail.com> wrote:
> 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.
>
> --
> 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*
More information about the gromacs.org_gmx-users
mailing list