[gmx-users] Free energy perturbation for larger molecule - big Coulomb solvation energies

[Krzysztof Makuch]

Hello!

I am trying to calculate solvation energy, using Free Energy Perturbation
for lutein in water and bilayer. Simulation runs smoothly, but Coulomb part
of calculated energy is extremly high for hydrophobic molecule
(C40H54(OH)2): above 40kcal/mol // 180 kJ/mol both in water and bilayer. I
am working with topology, with OPLS-AA charges.

To check possibility if the charges are somehow wrong I've parametrized
molecule with automatic server for charmm - cgenff. Parametrization of the
molecule in charmm clearly isn't perfect, but since what I want is only
rough estimation of possibility of wrong OPLS-AA charges it should be
sufficient. However, I am still getting even higher value of 300kJ/mol.
I've also checked methanol under same conditions in mdp file and obtained
results were reasonable.

Now, I suspect I've messed something in mdp FEP settings,could you please
look at it and give me a hint, what may be wrong with these?

Best,
Krzysiek

; Run control
integrator   = sd   ; Langevin dynamics
tinit= 0
dt   = 0.002
nsteps   = 175   ; 3.5 ns
nstcomm  = 100
; Output control
nstxout  = 500
nstvout  = 500
nstfout  = 0
nstlog   = 500
nstenergy= 500
nstxout-compressed   = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme= verlet
nstlist  = 20
ns_type  = grid
pbc  = xyz
rlist= 1.2
; Electrostatics
coulombtype  = PME
rcoulomb = 1.2
; van der Waals
vdwtype  = cutoff
vdw-modifier = potential-switch
rvdw-switch  = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing   = 0.1
; EWALD/PME/PPPM parameters
pme_order= 5
ewald_rtol   = 1e-06
epsilon_surface  = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps  = System
tau_t= 0.6
ref_t= 310
; Pressure coupling is on for NPT
Pcoupl   = Parrinello-Rahman
pcoupltype = isotropic
tau_p= 1.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_state012345
678910   11   12   13
14   15   16   17   18   19   20
vdw_lambdas  = 000000000
000000000000
coul_lambdas = 00.050.10.150.20.25
0.30.350.40.450.50.550.60.650.70.75
0.80.850.90.951
; We are not transforming any bonded or restrained interactions
bonded_lambdas   = 000000000
000000000000
restraint_lambdas= 000000000
000000000000
; Masses are not changing (particle identities are the same at lambda = 0
and lambda = 1)
mass_lambdas = 000000000
000000000000
; Not doing simulated temperting here
temperature_lambdas  = 000000000
000000000000
; Options for the decoupling
sc-alpha = 0.5
sc-coul  = no   ; linear interpolation of Coulomb (none
in this case)
sc-power = 1
sc-sigma = 0.3
couple-moltype   = lutein  ; name of moleculetype to decouple
couple-lambda0   = vdw-q  ; all
couple-lambda1   = vdw ; vdw
couple-intramol  = yes
nstdhdl  = 100
; No velocities during EM
gen_vel  = no
; 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  = 4

-- 
Jagiellonian University
Department of Computational Biophysics and Bioinformatics
tel.1: (12) 664 61 49
tel.2: (48) 664 086 049
-- 
Gromacs Users mailing list

* Please search the archive at 

Re: [gmx-users] antigen-antibody insilico MD using GROMACS

[Krzysztof Makuch]

Hi,
If you don't know where exactly binding occurs you can't really do this
using only gromacs. I'd even suggest that until you have to, for antibodies
it might be better to stick mostly to docking software.
1 - you have to have antibody and antigen parametrized using the same
forcefield (FF).
2 - check stability in selected FF. Some make proteins more stable, other
less. It is important to pay attention if your protein really have to be
stable - for example if you take only part of the antibody it may not be
stable either under physiological or in silico conditions. That's why you
may have to use some restrains. Or you may not, you have to decide.
3 - you have to dock antigen in antibody using docking software (there are
many...). Above 50% of results will be ridiculous, often the best energy
still means stupid docking result.
In many cases you may stop here, just use multiple docking engines, compare
and you have supplementary materials for your molecular biology paper. If
you are really stubborn and you have ligand parametrized:
4 - restrain both molecules, equilibrate
5 - turn off restrains and check if your complex is stable. For antibodies
there is quite big chance that it won't be because:
- Ab. is dynamic structure which may change shape after binding and docking
can't simulate this
- Water competes for H-bonds and nobody can guarantee that your complex is
stable enough
- you may have wrong var region.
- the position of antigen may be slightly wrong and in result - break off
7. If everything is ok you can begin energy analysis - umbrella, FEP
In general I advice against being stubborn until this is your primary
project. Antibodies are nasty, little guys.
Best,
KM

2018-02-23 15:24 GMT+01:00 Deep kumar :

> Hi Gromacs users,
>
> I am studying antigen-antibody interaction at protein level. I have a
> protein sequence (no crystal structure, length 180 residues) of the
> antigen, and have predicted the secondary structure of it (and modeled). By
> performing conserved domain search using inferno and NCBI I found out a
> domain in the antigen - PCC1 (residue 90 to 160). I am interested to know
> how the antigen interacts with the antibody (2000 residues) *in silico*.
> The ultimate aim is to find out "how the antigen interacts with the
> antibody, and predict the possible domains used by the antigen to interact
> with the antibody". I would really appreciate your valuable suggestions on
> how this can be done using GROMACS MD. Please let me know if I need to
> provide any further information. Thanks for your time.
>
> Regards,
> DK
> --
> 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-requ...@gromacs.org.
>



-- 
Jagiellonian University
Department of Computational Biophysics and Bioinformatics
tel.1: (12) 664 61 49
tel.2: (48) 664 086 049
-- 
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-requ...@gromacs.org.

[gmx-users] Solvation energy estimation gmx sasa

[Krzysztof Makuch]

Hi,
I'd like to ask about solvation energy given in gmx sasa (gromacs 2016)
with option:
-odg
What algorithm is used here, or rather - what values are given for area of
particular atom types? I guess there is publication behind this, but the
only mentioned in gmx sasa focus on solvation area. Is it the paper written
by Eisenberg & McLchlan in 1985?
I know this is only rough estimation, but still I'd love if it could be
mentioned in program description.
Best,
Krzysiek Makuch

-- 
Jagiellonian University
Department of Computational Biophysics and Bioinformatics
tel.1: (12) 664 61 49
tel.2: (48) 664 086 049
-- 
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-requ...@gromacs.org.

Re: [gmx-users] Do i need to put POSITION RESTRAINT DURING EQUILIBRATION STAGE ( NVT ) if i am preparing an amorphous sample?

[Krzysztof Makuch]

Hi,
That's always your call. You perform equilibration to swap between MM
(energy minimization) and MD (kinetic energy). In the other word you slowly
start MD and prevent unwanted and unrealistic rearrangement in your system.
If initial positions are important - you restrain the interesting
molecules. If you equilibrate for example lipid bilayer - there is no
reason to do that. Often you can also totally skip equilibration and just
cut off a few initial ns. Sometimes you make nvt, npt and still cut off
beginning and sometimes you perform MD to find optimal system, in which
case the whole simulation is in fact equilibration.
You know your system, you have to consider if restraining movement is what
you need.
Best,
KM

2018-02-09 13:28 GMT+01:00 sanjeet kumar singh ch16d012 <
ch16d...@smail.iitm.ac.in>:

> Hi list,
>
> I am preparing an amorphous sample using GROMACS but i am in doubt that
> during the equilibration stage ( NVT & NPT ) do i need to put position
> restraint on my polymer as there are no solvent in my system and if i have
> to use position restraint then why i should do that?
>
> THANKS,
> SK
> --
> 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-requ...@gromacs.org.
>



-- 
Jagiellonian University
Department of Computational Biophysics and Bioinformatics
tel.1: (12) 664 61 49
tel.2: (48) 664 086 049
-- 
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-requ...@gromacs.org.

Re: [gmx-users] Free energy perturbation - table extension and blowing system

[Krzysztof Makuch]

Ok, I've found what's wrong here. Intramolecular bonded interactions in
large molecule should be changed differently:
couple-intramol  = yes
Maybe somebody will find this useful.

2017-12-18 19:31 GMT+01:00 Krzysztof Makuch <krzysztof.mak...@gmail.com>:

> Hi,
> I want to calculate 'solvation' energy of lutein in bilayer and water by
> Free Energy Perturbation. When I'm running a normal MD, or some kind of
> steering dynamics,my system works just fine. After addition of FEP
> parameters, em runs, but instantly  returns series of warnings:
>
> WARNING: Listed nonbonded interaction between particles 6461 and 6482
> at distance 2.545 which is larger than the table limit 2.200 nm.
>
> This is likely either a 1,4 interaction, or a listed interaction inside
> a smaller molecule you are decoupling during a free energy calculation.
> Since interactions at distances beyond the table cannot be computed,
> they are skipped until they are inside the table limit again. You will
> only see this message once, even if it occurs for several interactions.
>
> IMPORTANT: This should not happen in a stable simulation, so there is
> probably something wrong with your system. Only change the table-extension
> distance in the mdp file if you are really sure that is the reason.
>
> The same thing repeats for several pairs of atoms, only in one
> transmembrane molecule. As I said, everything is fine without FEP
> parameters. For this instance I'm using the following parameters adapted
> stright from the Bevanlab tutorial:
>
> ; Run control
> integrator   = steep
> nsteps   = 5000
> ; EM criteria and other stuff
> emtol= 100
> emstep   = 0.01
> niter= 20
> nbfgscorr= 10
> ; Output control
> nstlog   = 1
> nstenergy= 1
> ; Neighborsearching and short-range nonbonded interactions
> cutoff-scheme= verlet
> nstlist  = 1
> ns_type  = grid
> pbc  = xyz
> rlist= 1.2
> ; Electrostatics
> coulombtype  = PME
> rcoulomb = 1.2
> ; van der Waals
> vdwtype  = cutoff
> vdw-modifier = potential-switch
> rvdw-switch  = 1.0
> rvdw = 1.2
> ; 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
> ; Temperature and pressure coupling are off during EM
> tcoupl   = no
> pcoupl   = no
> ; Free energy control stuff
> free_energy  = yes
> init_lambda_state= 0
> delta_lambda = 0
> calc_lambda_neighbors= 1
> ; Vectors of lambda specified here
> ; Each combination is an index that is retrieved from init_lambda_state
> for each simulation
> ; init_lambda_state012345678
> 910   11   12   13   14   15   16   17   18   19   20
> vdw_lambdas  = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
> 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
> coul_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
> ; 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
> sc-power = 1.0
> sc-sigma = 0.3
> couple-moltype   = lutein
> couple-lambda0   = vdw
> couple-lambda1   = none
> couple-intramol  = no
> nstdhdl  = 10
> ; No velocities during EM
> gen_vel  = no
> ; options for bonds
> constraints  = h-bonds
> ; Type of constraint algor

[gmx-users] Free energy perturbation - table extension and blowing system

[Krzysztof Makuch]

Hi,
I want to calculate 'solvation' energy of lutein in bilayer and water by
Free Energy Perturbation. When I'm running a normal MD, or some kind of
steering dynamics,my system works just fine. After addition of FEP
parameters, em runs, but instantly  returns series of warnings:

WARNING: Listed nonbonded interaction between particles 6461 and 6482
at distance 2.545 which is larger than the table limit 2.200 nm.

This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.

IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.

The same thing repeats for several pairs of atoms, only in one
transmembrane molecule. As I said, everything is fine without FEP
parameters. For this instance I'm using the following parameters adapted
stright from the Bevanlab tutorial:

; Run control
integrator   = steep
nsteps   = 5000
; EM criteria and other stuff
emtol= 100
emstep   = 0.01
niter= 20
nbfgscorr= 10
; Output control
nstlog   = 1
nstenergy= 1
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme= verlet
nstlist  = 1
ns_type  = grid
pbc  = xyz
rlist= 1.2
; Electrostatics
coulombtype  = PME
rcoulomb = 1.2
; van der Waals
vdwtype  = cutoff
vdw-modifier = potential-switch
rvdw-switch  = 1.0
rvdw = 1.2
; 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
; Temperature and pressure coupling are off during EM
tcoupl   = no
pcoupl   = no
; Free energy control stuff
free_energy  = yes
init_lambda_state= 0
delta_lambda = 0
calc_lambda_neighbors= 1
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for
each simulation
; init_lambda_state012345678
910   11   12   13   14   15   16   17   18   19   20
vdw_lambdas  = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_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
; 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
sc-power = 1.0
sc-sigma = 0.3
couple-moltype   = lutein
couple-lambda0   = vdw
couple-lambda1   = none
couple-intramol  = no
nstdhdl  = 10
; No velocities during EM
gen_vel  = no
; options for bonds
constraints  = h-bonds
; 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


What I've tried by now and didn't work:
- I've tried to increase table extension, but it didn't do any good.
- I've tried to change the charge groups, which are now 1 atom = 1 group
for this molecule
- I was trying to solve this with gromacs 4, 5 and 16'. Warnings were
similar.
If I try to run further FEP steps the system with FEP parameters blows up.
MD without FEP simulates doesn't cause any problems.
My guess is I'm missing something inside FEP parameters. I'm using GROMACS
for several years by now and for the first time I can't solve the problem
by myself.