Re: [gmx-users] Coarse-grained Protein-ligand simulations

[Justin Lemkul]

On 3/29/19 9:17 AM, Mac Kevin Braza wrote:

Thank you Professor Lemkul,

But would you suggest on how can I coarse-grained the ligand I am using? I
have been searching resources online but they do not work in our part.


I don't work with CG simulations, so I'm not much help. I would think 
that a CG parametrization of a ligand would remove all the detail you'd 
normally want to see in terms of ligand-protein interactions.


-Justin


I hope you can help us. Thank you Prof. Lemkul!

Best regards,
Mac Kevin E. Braza

On Fri, Mar 29, 2019, 8:59 PM Justin Lemkul  wrote:



On 3/29/19 3:32 AM, Mac Kevin Braza wrote:

Hello everyone,

I am simulating a coarse-grained model of a membrane protein (GPCR) in
lipid bilayer and an all-atom ligand octopamine. I build the protein,
solutes, and membrane in the web server CHARMM-GUI. While, I added the
ligand to the protein complex manually using the same coordinates of the
coarse-grained protein model.

I used the GROMACS input files from the output of CHARMM-GUI to simulate
the system. I include the LIGAND.ITP (from the PRODRG Server) to the
system.top and added the atom indexes in the index.ndx file.

Don't do this. An atomistic representation of a ligand and a CG
representation of everything else is incompatible. Mixing and matching
force fields is never a good idea. Moreover, PRODRG produces topologies
that are known to be unsuitable for MD simulations.


However, when I proceed with the second part of equilibration, the
following errors occurred.

*Command line*:
gmx grompp -f step6.2_equilibration.mdp -o step6.2_equilibration.tpr

-c

step6.1_equilibration.gro -p system.top -n index.ndx

Setting the LD random seed to 1722366284
Generated 2391 of the 4656 non-bonded parameter combinations
Excluding 1 bonded neighbours molecule type 'PROA_P'
Excluding 1 bonded neighbours molecule type 'POPC'
Excluding 1 bonded neighbours molecule type 'W'
Excluding 1 bonded neighbours molecule type 'NA'
Excluding 1 bonded neighbours molecule type 'CL'
Excluding 3 bonded neighbours molecule type 'LIG'
Velocities were taken from a Maxwell distribution at 303.15 K
Removing all charge groups because cutoff-scheme=Verlet

---
Program gmx grompp, VERSION 5.1.4
Source code file: /home/gromacs-5.1.4/src/gromacs/gmxpreprocess/readir.c,
line: 2690

Fatal error:
20 atoms are not part of any of the T-Coupling groups
For more information and tips for troubleshooting, please check the

GROMACS

website at http://www.gromacs.org/Documentation/Errors
---

The 20 atoms described the ligand I placed inside the protein-membrane
complex. I want to know if where can this error originate and how can we
fix them?

This simply means you haven't specified the ligand anywhere in tc-grps.
But again, back up and reevaluate your approach, which is far more
problematic than this simple index group issue.

-Justin

--
==

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalem...@vt.edu | (540) 231-3129
http://www.thelemkullab.com

==

--
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.



--
==

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalem...@vt.edu | (540) 231-3129
http://www.thelemkullab.com

==

--
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] Energy Conservation at the Beginning of a Production Run

[Justin Lemkul]

On 3/29/19 12:23 PM, Kruse, Luke E.(MU-Student) wrote:

Thank you all for the help!

Introducing the h-bond constraints into the minimization as well, did reduce 
the magnitude of the initial jump to 10.1 kJ per mole in the first 0.0066 ps 
but did not remove the problem entirely. I am also using different 
coulombtype's in the minimization (Cut-off) and the production run (PME). Could 
this be causing the spike, and, is this bad practice? The minimization did not 
converge when I used coulombtype = PME.


Cutoff electrostatics are horribly inaccurate and in and of themselves 
are bad practice. By changing the way the energy of the system is 
calculated between minimization and MD, you should certainly expect a 
difference in the energy. The same coordinates will, by definition, 
yield different energies.


-Justin



From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
 on behalf of Justin Lemkul 

Sent: Thursday, March 28, 2019 8:16:34 PM
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Energy Conservation at the Beginning of a Production 
Run



On 3/28/19 6:02 PM, Kruse, Luke E.(MU-Student) wrote:

In the production run I have constraints on h-bonds and I do not have any 
constraints in the minimization.

This would account for a sudden change in energy. At step 0, your H
positions will immediately change.

-Justin



From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
 on behalf of David van der Spoel 

Sent: Thursday, March 28, 2019 4:50:38 PM
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Energy Conservation at the Beginning of a Production 
Run

Den 2019-03-28 kl. 21:19, skrev Kruse, Luke E.(MU-Student):

In the mdout.mdp file corresponding to this run I have the line

; Do not constrain the start configuration
continuation = no


Do I need to change this to yes?

Maybe, try a 100 step run where you save energies at each time step. Did
you change anything else between minimization and production, e.g.
constraint settings?

The plot of energy vs time looks like an under damped response to a set point 
change in control theory - slight overshooting after the first few steps, then 
settling down to an average value nicely. This is what is leading me to believe 
that there is a large, non conservative force near the beginning that becomes 
smaller, and conservative.

Time (ps)

   0.00  -1851396.981075
   0.10  -1851388.994856
   0.20  -1851390.001591
   0.30  -1851389.436699
   0.40  -1851389.760231
   0.50  -1851389.726310
   0.60  -1851389.808428
   0.70  -1851389.907829
   0.80  -1851389.840644
   0.90  -1851389.967373
   1.00  -1851390.061502
   1.10  -1851390.057500
   1.20  -1851390.086172
   1.30  -1851390.042747
   1.40  -1851390.060210
   1.50  -1851389.998599
   1.60  -1851390.120319
   1.70  -1851390.141846
   1.80  -1851390.131003
   1.90  -1851390.170220
   2.00  -1851390.141617
   2.10  -1851390.135779
   2.20  -1851390.255565
   2.30  -1851390.153947
   2.40  -1851390.238681
   2.50  -1851390.144512
   2.60  -1851390.255200
   2.633000  -1851390.178012



From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
 on behalf of David van der Spoel 

Sent: Thursday, March 28, 2019 3:10:05 PM
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Energy Conservation at the Beginning of a Production 
Run

Den 2019-03-28 kl. 20:53, skrev Kruse, Luke E.(MU-Student):

Hello gromacs users,


I am trying to simulate a peptide amphiphile with the CHARMM27 force field. To 
do this I have had to specify an additional bond type and bond, angle, 
dihedral, etc. parameters in the .itp files of the force field. Then, to check 
if I had done this correctly, I minimized the system, and ran a production run 
to generate an NVE ensemble, so that I could make sure the system was 
conserving energy appropriately. After looking at the energy vs time plot 
(produces with the gmx energy command), however, the system jumps from an 
initial energy, up ~20 kJ per mole and then conserves energy for the most part 
(slight drifting but seems tolerable relative to the initial discontinuity). Is 
this discontinuity a normal happenstance or a result of bad minimization?


Is that at step 0 in the simulation?

Please check options like
; Do not constrain the start configuration
continuation = no



Thank you!

Luke


--
David van der Spoel, Ph.D., Professor of Biology
Head of Department, Cell & Molecular Biology, Uppsala University.
Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
http://www.icm.uu.se

Institutionen för cell- och molekylärbiologi - Uppsala 
universitet
www.icm.uu.se
Sanna Koskiniemi får ERC Starting Grant 2018-07-27 Ettemas 

Re: [gmx-users] Energy Conservation at the Beginning of a Production Run

[Kruse, Luke E.(MU-Student)]

Thank you all for the help!

Introducing the h-bond constraints into the minimization as well, did reduce 
the magnitude of the initial jump to 10.1 kJ per mole in the first 0.0066 ps 
but did not remove the problem entirely. I am also using different 
coulombtype's in the minimization (Cut-off) and the production run (PME). Could 
this be causing the spike, and, is this bad practice? The minimization did not 
converge when I used coulombtype = PME.


From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
 on behalf of Justin Lemkul 

Sent: Thursday, March 28, 2019 8:16:34 PM
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Energy Conservation at the Beginning of a Production 
Run



On 3/28/19 6:02 PM, Kruse, Luke E.(MU-Student) wrote:
> In the production run I have constraints on h-bonds and I do not have any 
> constraints in the minimization.

This would account for a sudden change in energy. At step 0, your H
positions will immediately change.

-Justin

> 
> From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
>  on behalf of David van 
> der Spoel 
> Sent: Thursday, March 28, 2019 4:50:38 PM
> To: gmx-us...@gromacs.org
> Subject: Re: [gmx-users] Energy Conservation at the Beginning of a Production 
> Run
>
> Den 2019-03-28 kl. 21:19, skrev Kruse, Luke E.(MU-Student):
>> In the mdout.mdp file corresponding to this run I have the line
>>
>> ; Do not constrain the start configuration
>> continuation = no
>>
>>
>> Do I need to change this to yes?
> Maybe, try a 100 step run where you save energies at each time step. Did
> you change anything else between minimization and production, e.g.
> constraint settings?
>>
>> The plot of energy vs time looks like an under damped response to a set 
>> point change in control theory - slight overshooting after the first few 
>> steps, then settling down to an average value nicely. This is what is 
>> leading me to believe that there is a large, non conservative force near the 
>> beginning that becomes smaller, and conservative.
>>
>> Time (ps)
>>
>>   0.00  -1851396.981075
>>   0.10  -1851388.994856
>>   0.20  -1851390.001591
>>   0.30  -1851389.436699
>>   0.40  -1851389.760231
>>   0.50  -1851389.726310
>>   0.60  -1851389.808428
>>   0.70  -1851389.907829
>>   0.80  -1851389.840644
>>   0.90  -1851389.967373
>>   1.00  -1851390.061502
>>   1.10  -1851390.057500
>>   1.20  -1851390.086172
>>   1.30  -1851390.042747
>>   1.40  -1851390.060210
>>   1.50  -1851389.998599
>>   1.60  -1851390.120319
>>   1.70  -1851390.141846
>>   1.80  -1851390.131003
>>   1.90  -1851390.170220
>>   2.00  -1851390.141617
>>   2.10  -1851390.135779
>>   2.20  -1851390.255565
>>   2.30  -1851390.153947
>>   2.40  -1851390.238681
>>   2.50  -1851390.144512
>>   2.60  -1851390.255200
>>   2.633000  -1851390.178012
>>
>>
>> 
>> From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
>>  on behalf of David van 
>> der Spoel 
>> Sent: Thursday, March 28, 2019 3:10:05 PM
>> To: gmx-us...@gromacs.org
>> Subject: Re: [gmx-users] Energy Conservation at the Beginning of a 
>> Production Run
>>
>> Den 2019-03-28 kl. 20:53, skrev Kruse, Luke E.(MU-Student):
>>> Hello gromacs users,
>>>
>>>
>>> I am trying to simulate a peptide amphiphile with the CHARMM27 force field. 
>>> To do this I have had to specify an additional bond type and bond, angle, 
>>> dihedral, etc. parameters in the .itp files of the force field. Then, to 
>>> check if I had done this correctly, I minimized the system, and ran a 
>>> production run to generate an NVE ensemble, so that I could make sure the 
>>> system was conserving energy appropriately. After looking at the energy vs 
>>> time plot (produces with the gmx energy command), however, the system jumps 
>>> from an initial energy, up ~20 kJ per mole and then conserves energy for 
>>> the most part (slight drifting but seems tolerable relative to the initial 
>>> discontinuity). Is this discontinuity a normal happenstance or a result of 
>>> bad minimization?
>>>
>> Is that at step 0 in the simulation?
>>
>> Please check options like
>> ; Do not constrain the start configuration
>> continuation = no
>>
>>
>>> Thank you!
>>>
>>> Luke
>>>
>>
>> --
>> David van der Spoel, Ph.D., Professor of Biology
>> Head of Department, Cell & Molecular Biology, Uppsala University.
>> Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
>> http://www.icm.uu.se
Institutionen för cell- och molekylärbiologi - Uppsala 
universitet
www.icm.uu.se
Sanna Koskiniemi får ERC Starting Grant 2018-07-27 Ettemas grupp klargör 
mitokondriellt ursprung i Nature journal 2018-05-25 Johanssons och Elfs grupp 
rapporterar om ny metod för 

Re: [gmx-users] PME constant shifts.

[Sergio Perez]

Great! So I just have to discount this term to obtain the positive energy
expected.
Thank you very much!
Sergio

On Fri, Mar 29, 2019 at 4:11 PM Berk Hess  wrote:

> I mean the q_tot^2 term.
> This attractative term goes with the square of the total charge, so has
> larger magnitude for AB than A+B, so the result in negative.
>
> Berk
>
> 
> From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <
> gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of Sergio
> Perez 
> Sent: Friday, March 29, 2019 3:15 PM
> To: gmx-us...@gromacs.org
> Subject: Re: [gmx-users] PME constant shifts.
>
> Hello,
>
> Thanks for the quick reply, I am not sure if you are referring to this term
> of the Ewald Sum:
>
> V_0 = ( f*beta ) * SUM(q_i^2) / sqrt(pi)
>
> Since the sum of the square of the charges of the system (AB) is the same
> as the sum of the square of the charges of the system (A) plus the same
> from system (B), this term should vanish.
>
> There is this an additional term for the Ewald Sum when the system is
> charged to compensate for the uniform charge as described by K. Smith in
> Elements of Molecular Dynamics section 6.4.6 (free book online):
>
> V= - (q_tot^2)/(8*eps_0*V*Beta)
>
> If this was not included by gromacs then this would explain why I get
> negative interaction energies.
>
> I know charged systems should be avoided (for example if the system is
> inhomogeneous, as grompp points out) but my purposes are FF development and
> finding bugs in my FF. Although the absolute energies are meaningless it
> would be nice if they could be intuitively correct.
>
> Best regards and thanks again!
>
> Sergio
>
>
> On Fri, Mar 29, 2019 at 1:45 PM Berk Hess  wrote:
>
> > Hi,
> >
> > The negative energy is due to the uniform background charge  that appears
> > when your system has a net charge, not due to shifting of pair
> > interactions. Computing the absolute energy of a periodic system with net
> > charge is meaningless.
> >
> > Cheers,
> >
> > Berk
> >
> > 
> > From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <
> > gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of Sergio
> > Perez 
> > Sent: Friday, March 29, 2019 12:33 PM
> > To: gromacs.org_gmx-users@maillist.sys.kth.se
> > Subject: [gmx-users] PME constant shifts.
> >
> > Dear gmx comunity,
> >
> > I am trying to calculate the electrostatic interaction of my system for
> > force field development. My system consists of 7 positively charged
> > particles interacting (system A) with another positively charged particle
> > (system B). I calculate the electrostatic interaction by doing reruns of
> > the trajectory in which some of them have been removed.
> >
> > E_int(electrostatic) =
> >
> >
> E_AB(coul-SR)+E_AB(coul-LR)-E_A(coul-SR)-E_A(coul-LR)-E_B(coul-SR)-E_B(coul-LR)
> >
> > And to my surprise I get negative energies. I know the PME method has its
> > energy shifted:
> >
> > "With the Verlet cut-off scheme, the PME direct space potential is
> shifted
> > by a constant such that
> > the potential is zero at the cut-off. This shift is small and since the
> net
> > system charge is close to
> > zero, the total shift is very small, unlike in the case of the
> > Lennard-Jones potential where all shifts
> > add up. We apply the shift anyhow, such that the potential is the exact
> > integral of the force."
> >
> > Since for the system with just one charge I get a possitive
> E_1q(coul-LR),
> > I imagine the LR-coul term is shifted (eventhough I have to use the group
> > cut-off scheme and not the verlet). The problem is why don't all these
> > shifts cancel? Is there a way to not add the shifts?
> >
> > I leave bellow the .mdp parameters. Note that for the VDW part I need to
> > use user-tables and therefore I need to use the group scheme.
> >
> > Thank you very much!
> > Sergio Pérez-Conesa
> >
> >
> > ; Integrator stuff
> > integrator = md-vv   ; steep or md
> > dt = 0.001
> > nsteps = 3 ; -1 no maximum
> > init-step = 30305000 ;
> > ld-seed = -2018962629
> > continuation = yes
> > emtol   = 10.0
> > emstep  = 0.01
> >
> > ; neighbour search
> > cutoff-scheme = group
> > nst-list = 10
> > verlet-buffer-tolerance = 0.005
> > ns-type = grid
> > ; Pbc
> > pbc = xyz ; xyz or xy
> >
> > ; Vdw
> > rvdw = 1.4
> > rlist = 1.4
> > vdwtype = user  ; PME or User to look for table /user
> > DispCorr = No   ; corrections to energy and pressure or No
> >
> > ;Electrostatic
> > coulombtype = PME   ; User if look for table
> > rcoulomb = 1.4
> > pme-order =  4
> > fourierspacing = 0.2
> > ewald-geometry = 3d ; 3d or 3dc
> >
> > ;IMD-group=
> >
> > ; T and P
> > tcoupl =no  ; none, nose-hoover, v-rescale
> > nh-chain-length = 10
> > nsttcouple = 10 ; 10 if not used
> > pcoupl = no ;
> > pcoupltype = semiisotropic ;
> > compressibility = 4.5e-5 4.5e-5  ; isothermal compressibility of
> > water, bar^-1
> > refcoord-scaling = com
> > tau-p = 

Re: [gmx-users] PME constant shifts.

[Berk Hess]

I mean the q_tot^2 term.
This attractative term goes with the square of the total charge, so has larger 
magnitude for AB than A+B, so the result in negative.

Berk


From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
 on behalf of Sergio Perez 

Sent: Friday, March 29, 2019 3:15 PM
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] PME constant shifts.

Hello,

Thanks for the quick reply, I am not sure if you are referring to this term
of the Ewald Sum:

V_0 = ( f*beta ) * SUM(q_i^2) / sqrt(pi)

Since the sum of the square of the charges of the system (AB) is the same
as the sum of the square of the charges of the system (A) plus the same
from system (B), this term should vanish.

There is this an additional term for the Ewald Sum when the system is
charged to compensate for the uniform charge as described by K. Smith in
Elements of Molecular Dynamics section 6.4.6 (free book online):

V= - (q_tot^2)/(8*eps_0*V*Beta)

If this was not included by gromacs then this would explain why I get
negative interaction energies.

I know charged systems should be avoided (for example if the system is
inhomogeneous, as grompp points out) but my purposes are FF development and
finding bugs in my FF. Although the absolute energies are meaningless it
would be nice if they could be intuitively correct.

Best regards and thanks again!

Sergio


On Fri, Mar 29, 2019 at 1:45 PM Berk Hess  wrote:

> Hi,
>
> The negative energy is due to the uniform background charge  that appears
> when your system has a net charge, not due to shifting of pair
> interactions. Computing the absolute energy of a periodic system with net
> charge is meaningless.
>
> Cheers,
>
> Berk
>
> 
> From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <
> gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of Sergio
> Perez 
> Sent: Friday, March 29, 2019 12:33 PM
> To: gromacs.org_gmx-users@maillist.sys.kth.se
> Subject: [gmx-users] PME constant shifts.
>
> Dear gmx comunity,
>
> I am trying to calculate the electrostatic interaction of my system for
> force field development. My system consists of 7 positively charged
> particles interacting (system A) with another positively charged particle
> (system B). I calculate the electrostatic interaction by doing reruns of
> the trajectory in which some of them have been removed.
>
> E_int(electrostatic) =
>
> E_AB(coul-SR)+E_AB(coul-LR)-E_A(coul-SR)-E_A(coul-LR)-E_B(coul-SR)-E_B(coul-LR)
>
> And to my surprise I get negative energies. I know the PME method has its
> energy shifted:
>
> "With the Verlet cut-off scheme, the PME direct space potential is shifted
> by a constant such that
> the potential is zero at the cut-off. This shift is small and since the net
> system charge is close to
> zero, the total shift is very small, unlike in the case of the
> Lennard-Jones potential where all shifts
> add up. We apply the shift anyhow, such that the potential is the exact
> integral of the force."
>
> Since for the system with just one charge I get a possitive E_1q(coul-LR),
> I imagine the LR-coul term is shifted (eventhough I have to use the group
> cut-off scheme and not the verlet). The problem is why don't all these
> shifts cancel? Is there a way to not add the shifts?
>
> I leave bellow the .mdp parameters. Note that for the VDW part I need to
> use user-tables and therefore I need to use the group scheme.
>
> Thank you very much!
> Sergio Pérez-Conesa
>
>
> ; Integrator stuff
> integrator = md-vv   ; steep or md
> dt = 0.001
> nsteps = 3 ; -1 no maximum
> init-step = 30305000 ;
> ld-seed = -2018962629
> continuation = yes
> emtol   = 10.0
> emstep  = 0.01
>
> ; neighbour search
> cutoff-scheme = group
> nst-list = 10
> verlet-buffer-tolerance = 0.005
> ns-type = grid
> ; Pbc
> pbc = xyz ; xyz or xy
>
> ; Vdw
> rvdw = 1.4
> rlist = 1.4
> vdwtype = user  ; PME or User to look for table /user
> DispCorr = No   ; corrections to energy and pressure or No
>
> ;Electrostatic
> coulombtype = PME   ; User if look for table
> rcoulomb = 1.4
> pme-order =  4
> fourierspacing = 0.2
> ewald-geometry = 3d ; 3d or 3dc
>
> ;IMD-group=
>
> ; T and P
> tcoupl =no  ; none, nose-hoover, v-rescale
> nh-chain-length = 10
> nsttcouple = 10 ; 10 if not used
> pcoupl = no ;
> pcoupltype = semiisotropic ;
> compressibility = 4.5e-5 4.5e-5  ; isothermal compressibility of
> water, bar^-1
> refcoord-scaling = com
> tau-p = 5.0  ; ps
> ref-p = 1.0 1.0 ; bar
>
> ; Constraints
> constraints = none  ; constraints  distances
> constraint-algorithm = LINCS ; or SHAKE
> lincs_iter  = 1
> lincs_order = 4
>
> ;Generate velocities
> gen_vel = no; assign velocities from Maxwell
> distribution
> gen_temp= 300   ; temperature for Maxwell distribution
> gen_seed= 180611; generate a random seed
>
> ; OUTPUT
> ; Output control
> 

Re: [gmx-users] PME constant shifts.

[Sergio Perez]

Hello,

Thanks for the quick reply, I am not sure if you are referring to this term
of the Ewald Sum:

V_0 = ( f*beta ) * SUM(q_i^2) / sqrt(pi)

Since the sum of the square of the charges of the system (AB) is the same
as the sum of the square of the charges of the system (A) plus the same
from system (B), this term should vanish.

There is this an additional term for the Ewald Sum when the system is
charged to compensate for the uniform charge as described by K. Smith in
Elements of Molecular Dynamics section 6.4.6 (free book online):

V= - (q_tot^2)/(8*eps_0*V*Beta)

If this was not included by gromacs then this would explain why I get
negative interaction energies.

I know charged systems should be avoided (for example if the system is
inhomogeneous, as grompp points out) but my purposes are FF development and
finding bugs in my FF. Although the absolute energies are meaningless it
would be nice if they could be intuitively correct.

Best regards and thanks again!

Sergio


On Fri, Mar 29, 2019 at 1:45 PM Berk Hess  wrote:

> Hi,
>
> The negative energy is due to the uniform background charge  that appears
> when your system has a net charge, not due to shifting of pair
> interactions. Computing the absolute energy of a periodic system with net
> charge is meaningless.
>
> Cheers,
>
> Berk
>
> 
> From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <
> gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of Sergio
> Perez 
> Sent: Friday, March 29, 2019 12:33 PM
> To: gromacs.org_gmx-users@maillist.sys.kth.se
> Subject: [gmx-users] PME constant shifts.
>
> Dear gmx comunity,
>
> I am trying to calculate the electrostatic interaction of my system for
> force field development. My system consists of 7 positively charged
> particles interacting (system A) with another positively charged particle
> (system B). I calculate the electrostatic interaction by doing reruns of
> the trajectory in which some of them have been removed.
>
> E_int(electrostatic) =
>
> E_AB(coul-SR)+E_AB(coul-LR)-E_A(coul-SR)-E_A(coul-LR)-E_B(coul-SR)-E_B(coul-LR)
>
> And to my surprise I get negative energies. I know the PME method has its
> energy shifted:
>
> "With the Verlet cut-off scheme, the PME direct space potential is shifted
> by a constant such that
> the potential is zero at the cut-off. This shift is small and since the net
> system charge is close to
> zero, the total shift is very small, unlike in the case of the
> Lennard-Jones potential where all shifts
> add up. We apply the shift anyhow, such that the potential is the exact
> integral of the force."
>
> Since for the system with just one charge I get a possitive E_1q(coul-LR),
> I imagine the LR-coul term is shifted (eventhough I have to use the group
> cut-off scheme and not the verlet). The problem is why don't all these
> shifts cancel? Is there a way to not add the shifts?
>
> I leave bellow the .mdp parameters. Note that for the VDW part I need to
> use user-tables and therefore I need to use the group scheme.
>
> Thank you very much!
> Sergio Pérez-Conesa
>
>
> ; Integrator stuff
> integrator = md-vv   ; steep or md
> dt = 0.001
> nsteps = 3 ; -1 no maximum
> init-step = 30305000 ;
> ld-seed = -2018962629
> continuation = yes
> emtol   = 10.0
> emstep  = 0.01
>
> ; neighbour search
> cutoff-scheme = group
> nst-list = 10
> verlet-buffer-tolerance = 0.005
> ns-type = grid
> ; Pbc
> pbc = xyz ; xyz or xy
>
> ; Vdw
> rvdw = 1.4
> rlist = 1.4
> vdwtype = user  ; PME or User to look for table /user
> DispCorr = No   ; corrections to energy and pressure or No
>
> ;Electrostatic
> coulombtype = PME   ; User if look for table
> rcoulomb = 1.4
> pme-order =  4
> fourierspacing = 0.2
> ewald-geometry = 3d ; 3d or 3dc
>
> ;IMD-group=
>
> ; T and P
> tcoupl =no  ; none, nose-hoover, v-rescale
> nh-chain-length = 10
> nsttcouple = 10 ; 10 if not used
> pcoupl = no ;
> pcoupltype = semiisotropic ;
> compressibility = 4.5e-5 4.5e-5  ; isothermal compressibility of
> water, bar^-1
> refcoord-scaling = com
> tau-p = 5.0  ; ps
> ref-p = 1.0 1.0 ; bar
>
> ; Constraints
> constraints = none  ; constraints  distances
> constraint-algorithm = LINCS ; or SHAKE
> lincs_iter  = 1
> lincs_order = 4
>
> ;Generate velocities
> gen_vel = no; assign velocities from Maxwell
> distribution
> gen_temp= 300   ; temperature for Maxwell distribution
> gen_seed= 180611; generate a random seed
>
> ; OUTPUT
> ; Output control
> nstenergy   = 1 ; save energies
> nstlog  = 10 ; update log file every
> ;nstxout-compressed  = 500  ; save compressed coordinates every 10.0 ps
> ;nstxout = 20
> ;nstvout = 20
> ;nstfout = 20
>
> ; COM removal . There are more options like removing from groups or the
> frequency
> ;comm-grps = System
> comm-mode = none; linear, angular (removes linear and 

Re: [gmx-users] Coarse-grained Protein-ligand simulations

[Mac Kevin Braza]

Thank you Professor Lemkul,

But would you suggest on how can I coarse-grained the ligand I am using? I
have been searching resources online but they do not work in our part.

I hope you can help us. Thank you Prof. Lemkul!

Best regards,
Mac Kevin E. Braza

On Fri, Mar 29, 2019, 8:59 PM Justin Lemkul  wrote:

>
>
> On 3/29/19 3:32 AM, Mac Kevin Braza wrote:
> > Hello everyone,
> >
> > I am simulating a coarse-grained model of a membrane protein (GPCR) in
> > lipid bilayer and an all-atom ligand octopamine. I build the protein,
> > solutes, and membrane in the web server CHARMM-GUI. While, I added the
> > ligand to the protein complex manually using the same coordinates of the
> > coarse-grained protein model.
> >
> > I used the GROMACS input files from the output of CHARMM-GUI to simulate
> > the system. I include the LIGAND.ITP (from the PRODRG Server) to the
> > system.top and added the atom indexes in the index.ndx file.
>
> Don't do this. An atomistic representation of a ligand and a CG
> representation of everything else is incompatible. Mixing and matching
> force fields is never a good idea. Moreover, PRODRG produces topologies
> that are known to be unsuitable for MD simulations.
>
> > However, when I proceed with the second part of equilibration, the
> > following errors occurred.
> >
> > *Command line*:
> >gmx grompp -f step6.2_equilibration.mdp -o step6.2_equilibration.tpr
> -c
> > step6.1_equilibration.gro -p system.top -n index.ndx
> >
> > Setting the LD random seed to 1722366284
> > Generated 2391 of the 4656 non-bonded parameter combinations
> > Excluding 1 bonded neighbours molecule type 'PROA_P'
> > Excluding 1 bonded neighbours molecule type 'POPC'
> > Excluding 1 bonded neighbours molecule type 'W'
> > Excluding 1 bonded neighbours molecule type 'NA'
> > Excluding 1 bonded neighbours molecule type 'CL'
> > Excluding 3 bonded neighbours molecule type 'LIG'
> > Velocities were taken from a Maxwell distribution at 303.15 K
> > Removing all charge groups because cutoff-scheme=Verlet
> >
> > ---
> > Program gmx grompp, VERSION 5.1.4
> > Source code file: /home/gromacs-5.1.4/src/gromacs/gmxpreprocess/readir.c,
> > line: 2690
> >
> > Fatal error:
> > 20 atoms are not part of any of the T-Coupling groups
> > For more information and tips for troubleshooting, please check the
> GROMACS
> > website at http://www.gromacs.org/Documentation/Errors
> > ---
> >
> > The 20 atoms described the ligand I placed inside the protein-membrane
> > complex. I want to know if where can this error originate and how can we
> > fix them?
>
> This simply means you haven't specified the ligand anywhere in tc-grps.
> But again, back up and reevaluate your approach, which is far more
> problematic than this simple index group issue.
>
> -Justin
>
> --
> ==
>
> Justin A. Lemkul, Ph.D.
> Assistant Professor
> Office: 301 Fralin Hall
> Lab: 303 Engel Hall
>
> Virginia Tech Department of Biochemistry
> 340 West Campus Dr.
> Blacksburg, VA 24061
>
> jalem...@vt.edu | (540) 231-3129
> http://www.thelemkullab.com
>
> ==
>
> --
> 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.
>
-- 
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] Help needed in using gmx select

[Budheswar Dehury]

Dear Gromacs Developers and Users,
I need a small help in GROMACS trajectory analysis. I want to define a 
selection i.e., 4Ang volume around two residues (from different chains X  and 
chain Y ) that make your binding site --- Within this selection I want count 
number of waters molecules along the trajectory (it's a lipid bilayer system 
with TIP3P and POPC).

How can I estimate the number of water molecules using gmx select and gmx dist 
or gmx mindist..through gromacs or VMD.

Looking forward to hear from you

Thanking You
Sincere Regards

Budheswar








-- 
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] Energy Conservation at the Beginning of a Production Run

[Justin Lemkul]

On 3/29/19 4:30 AM, Mark Abraham wrote:

Hi,

Off topic for this discussion, but what would people think if we changed
things around so that grompp would do that constraining, before writing the
tpr file?

This will mean that the configuration in the tpr does not always match the
configuration that was given to grompp, but instead matches the rest of the
simulation. Thus such an energy jump would not occur. And it may mean that
broken structures are found at grompp time rather than later at step 0 when
doing mdrun on the cluster.


I would think any problems that might arise would be minimal, such that 
the user can easily recover the original geometry by turning constraints 
off and generating a new .tpr file. This could also be a very good 
opportunity to notify the user if they are perhaps unwittingly 
introducing a discontinuity between unconstrained and constrained 
processes, e.g. if the constraint RMSD is larger than some (small) 
threshold, we can print a note to alert the user about potential 
discontinuities.


-Justin


Mark

On Fri., 29 Mar. 2019, 02:17 Justin Lemkul,  wrote:



On 3/28/19 6:02 PM, Kruse, Luke E.(MU-Student) wrote:

In the production run I have constraints on h-bonds and I do not have

any constraints in the minimization.

This would account for a sudden change in energy. At step 0, your H
positions will immediately change.

-Justin



From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <

gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of David van
der Spoel 

Sent: Thursday, March 28, 2019 4:50:38 PM
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Energy Conservation at the Beginning of a

Production Run

Den 2019-03-28 kl. 21:19, skrev Kruse, Luke E.(MU-Student):

In the mdout.mdp file corresponding to this run I have the line

; Do not constrain the start configuration
continuation = no


Do I need to change this to yes?

Maybe, try a 100 step run where you save energies at each time step. Did
you change anything else between minimization and production, e.g.
constraint settings?

The plot of energy vs time looks like an under damped response to a set

point change in control theory - slight overshooting after the first few
steps, then settling down to an average value nicely. This is what is
leading me to believe that there is a large, non conservative force near
the beginning that becomes smaller, and conservative.

Time (ps)

   0.00  -1851396.981075
   0.10  -1851388.994856
   0.20  -1851390.001591
   0.30  -1851389.436699
   0.40  -1851389.760231
   0.50  -1851389.726310
   0.60  -1851389.808428
   0.70  -1851389.907829
   0.80  -1851389.840644
   0.90  -1851389.967373
   1.00  -1851390.061502
   1.10  -1851390.057500
   1.20  -1851390.086172
   1.30  -1851390.042747
   1.40  -1851390.060210
   1.50  -1851389.998599
   1.60  -1851390.120319
   1.70  -1851390.141846
   1.80  -1851390.131003
   1.90  -1851390.170220
   2.00  -1851390.141617
   2.10  -1851390.135779
   2.20  -1851390.255565
   2.30  -1851390.153947
   2.40  -1851390.238681
   2.50  -1851390.144512
   2.60  -1851390.255200
   2.633000  -1851390.178012



From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <

gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of David van
der Spoel 

Sent: Thursday, March 28, 2019 3:10:05 PM
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Energy Conservation at the Beginning of a

Production Run

Den 2019-03-28 kl. 20:53, skrev Kruse, Luke E.(MU-Student):

Hello gromacs users,


I am trying to simulate a peptide amphiphile with the CHARMM27 force

field. To do this I have had to specify an additional bond type and bond,
angle, dihedral, etc. parameters in the .itp files of the force field.
Then, to check if I had done this correctly, I minimized the system, and
ran a production run to generate an NVE ensemble, so that I could make sure
the system was conserving energy appropriately. After looking at the energy
vs time plot (produces with the gmx energy command), however, the system
jumps from an initial energy, up ~20 kJ per mole and then conserves energy
for the most part (slight drifting but seems tolerable relative to the
initial discontinuity). Is this discontinuity a normal happenstance or a
result of bad minimization?

Is that at step 0 in the simulation?

Please check options like
; Do not constrain the start configuration
continuation = no



Thank you!

Luke


--
David van der Spoel, Ph.D., Professor of Biology
Head of Department, Cell & Molecular Biology, Uppsala University.
Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
http://www.icm.uu.se
--
Gromacs Users mailing list

* Please search the archive at

Re: [gmx-users] Coarse-grained Protein-ligand simulations

[Justin Lemkul]

On 3/29/19 3:32 AM, Mac Kevin Braza wrote:

Hello everyone,

I am simulating a coarse-grained model of a membrane protein (GPCR) in
lipid bilayer and an all-atom ligand octopamine. I build the protein,
solutes, and membrane in the web server CHARMM-GUI. While, I added the
ligand to the protein complex manually using the same coordinates of the
coarse-grained protein model.

I used the GROMACS input files from the output of CHARMM-GUI to simulate
the system. I include the LIGAND.ITP (from the PRODRG Server) to the
system.top and added the atom indexes in the index.ndx file.


Don't do this. An atomistic representation of a ligand and a CG 
representation of everything else is incompatible. Mixing and matching 
force fields is never a good idea. Moreover, PRODRG produces topologies 
that are known to be unsuitable for MD simulations.



However, when I proceed with the second part of equilibration, the
following errors occurred.

*Command line*:
   gmx grompp -f step6.2_equilibration.mdp -o step6.2_equilibration.tpr -c
step6.1_equilibration.gro -p system.top -n index.ndx

Setting the LD random seed to 1722366284
Generated 2391 of the 4656 non-bonded parameter combinations
Excluding 1 bonded neighbours molecule type 'PROA_P'
Excluding 1 bonded neighbours molecule type 'POPC'
Excluding 1 bonded neighbours molecule type 'W'
Excluding 1 bonded neighbours molecule type 'NA'
Excluding 1 bonded neighbours molecule type 'CL'
Excluding 3 bonded neighbours molecule type 'LIG'
Velocities were taken from a Maxwell distribution at 303.15 K
Removing all charge groups because cutoff-scheme=Verlet

---
Program gmx grompp, VERSION 5.1.4
Source code file: /home/gromacs-5.1.4/src/gromacs/gmxpreprocess/readir.c,
line: 2690

Fatal error:
20 atoms are not part of any of the T-Coupling groups
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---

The 20 atoms described the ligand I placed inside the protein-membrane
complex. I want to know if where can this error originate and how can we
fix them?


This simply means you haven't specified the ligand anywhere in tc-grps. 
But again, back up and reevaluate your approach, which is far more 
problematic than this simple index group issue.


-Justin

--
==

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalem...@vt.edu | (540) 231-3129
http://www.thelemkullab.com

==

--
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] PME constant shifts. (Sergio Perez)

[Groenhof, Gerrit]

Dear Sergio,


If you're using PME, you need to also include the reciprocal contribution.


Gerrit



Dear gmx comunity,

I am trying to calculate the electrostatic interaction of my system for
force field development. My system consists of 7 positively charged
particles interacting (system A) with another positively charged particle
(system B). I calculate the electrostatic interaction by doing reruns of
the trajectory in which some of them have been removed.

E_int(electrostatic) =
E_AB(coul-SR)+E_AB(coul-LR)-E_A(coul-SR)-E_A(coul-LR)-E_B(coul-SR)-E_B(coul-LR)

And to my surprise I get negative energies. I know the PME method has its
energy shifted:

"With the Verlet cut-off scheme, the PME direct space potential is shifted
by a constant such that
the potential is zero at the cut-off. This shift is small and since the net
system charge is close to
zero, the total shift is very small, unlike in the case of the
Lennard-Jones potential where all shifts
add up. We apply the shift anyhow, such that the potential is the exact
integral of the force."

Since for the system with just one charge I get a possitive E_1q(coul-LR),
I imagine the LR-coul term is shifted (eventhough I have to use the group
cut-off scheme and not the verlet). The problem is why don't all these
shifts cancel? Is there a way to not add the shifts?

I leave bellow the .mdp parameters. Note that for the VDW part I need to
use user-tables and therefore I need to use the group scheme.

Thank you very much!
Sergio P?rez-Conesa


; neighbour search
cutoff-scheme = group
nst-list = 10
verlet-buffer-tolerance = 0.005
ns-type = grid
; Pbc
pbc = xyz ; xyz or xy

; Vdw
rvdw = 1.4
rlist = 1.4
vdwtype = user  ; PME or User to look for table /user
DispCorr = No   ; corrections to energy and pressure or No

;Electrostatic
coulombtype = PME   ; User if look for table
rcoulomb = 1.4
pme-order =  4
fourierspacing = 0.2
ewald-geometry = 3d ; 3d or 3dc

;IMD-
-- 
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] CG lipid diffusion coefficient

[Yasser Almeida Hernández]

Hi all,

I ran some CG simulations of a protein inserted in a lipid bilayer, with 
several lipids. I used the msd tool to compute the lipid (using the PO4) 
diffusion coefficient. For one of the lipid I got D = 0.6 x 10^-5 
cm^2/s, which is 600 um^2/s. According to the Martini website "To 
compare the diffusion coefficient obtained from a Martini simulation to 
a measured one, a conversion factor of about 4 has to be applied to 
account for the faster diffusion at the CG level due to the smoothened 
free energy landscape (note, the conversion factor can vary 
significantly depending on the molecule in question)." (J Phys Chem B. 
2009 Feb 5; 113(5): 1501–1510.)


The value I get is 2 orders of magnitude higher of a comparable value of 
6 um^2/s.


Would this make sense?

Regards

Yasser

--
Yasser Almeida Hernández
PhD student
Institute of Biochemistry and Molecular Biology
Department of Chemistry
University of Hamburg
Martin-Luther-King-Platz 6
20146 Hamburg
Germany
+49 40 42838 2845
yasser.almeida.hernan...@chemie.uni-hamburg.de
office: Grindelallee 117, room 250c

--
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] PME constant shifts.

[Berk Hess]

Hi,

The negative energy is due to the uniform background charge  that appears when 
your system has a net charge, not due to shifting of pair interactions. 
Computing the absolute energy of a periodic system with net charge is 
meaningless.

Cheers,

Berk


From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
 on behalf of Sergio Perez 

Sent: Friday, March 29, 2019 12:33 PM
To: gromacs.org_gmx-users@maillist.sys.kth.se
Subject: [gmx-users] PME constant shifts.

Dear gmx comunity,

I am trying to calculate the electrostatic interaction of my system for
force field development. My system consists of 7 positively charged
particles interacting (system A) with another positively charged particle
(system B). I calculate the electrostatic interaction by doing reruns of
the trajectory in which some of them have been removed.

E_int(electrostatic) =
E_AB(coul-SR)+E_AB(coul-LR)-E_A(coul-SR)-E_A(coul-LR)-E_B(coul-SR)-E_B(coul-LR)

And to my surprise I get negative energies. I know the PME method has its
energy shifted:

"With the Verlet cut-off scheme, the PME direct space potential is shifted
by a constant such that
the potential is zero at the cut-off. This shift is small and since the net
system charge is close to
zero, the total shift is very small, unlike in the case of the
Lennard-Jones potential where all shifts
add up. We apply the shift anyhow, such that the potential is the exact
integral of the force."

Since for the system with just one charge I get a possitive E_1q(coul-LR),
I imagine the LR-coul term is shifted (eventhough I have to use the group
cut-off scheme and not the verlet). The problem is why don't all these
shifts cancel? Is there a way to not add the shifts?

I leave bellow the .mdp parameters. Note that for the VDW part I need to
use user-tables and therefore I need to use the group scheme.

Thank you very much!
Sergio Pérez-Conesa


; Integrator stuff
integrator = md-vv   ; steep or md
dt = 0.001
nsteps = 3 ; -1 no maximum
init-step = 30305000 ;
ld-seed = -2018962629
continuation = yes
emtol   = 10.0
emstep  = 0.01

; neighbour search
cutoff-scheme = group
nst-list = 10
verlet-buffer-tolerance = 0.005
ns-type = grid
; Pbc
pbc = xyz ; xyz or xy

; Vdw
rvdw = 1.4
rlist = 1.4
vdwtype = user  ; PME or User to look for table /user
DispCorr = No   ; corrections to energy and pressure or No

;Electrostatic
coulombtype = PME   ; User if look for table
rcoulomb = 1.4
pme-order =  4
fourierspacing = 0.2
ewald-geometry = 3d ; 3d or 3dc

;IMD-group=

; T and P
tcoupl =no  ; none, nose-hoover, v-rescale
nh-chain-length = 10
nsttcouple = 10 ; 10 if not used
pcoupl = no ;
pcoupltype = semiisotropic ;
compressibility = 4.5e-5 4.5e-5  ; isothermal compressibility of
water, bar^-1
refcoord-scaling = com
tau-p = 5.0  ; ps
ref-p = 1.0 1.0 ; bar

; Constraints
constraints = none  ; constraints  distances
constraint-algorithm = LINCS ; or SHAKE
lincs_iter  = 1
lincs_order = 4

;Generate velocities
gen_vel = no; assign velocities from Maxwell
distribution
gen_temp= 300   ; temperature for Maxwell distribution
gen_seed= 180611; generate a random seed

; OUTPUT
; Output control
nstenergy   = 1 ; save energies
nstlog  = 10 ; update log file every
;nstxout-compressed  = 500  ; save compressed coordinates every 10.0 ps
;nstxout = 20
;nstvout = 20
;nstfout = 20

; COM removal . There are more options like removing from groups or the
frequency
;comm-grps = System
comm-mode = none; linear, angular (removes linear and angular motion)
or none
nstcomm = 0  ; frequency of removal of com motion
freezegrps =
energygrp-excl =
;freezedim = Y Y Y
nstcalcenergy =
--
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.
-- 
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] WG: WG: Issue with CUDA and gromacs

[Szilárd Páll]

Hi,

The standard output of the first set of runs is also something I was
interested in, but I've found the equivalent in the
complex/TESTDIR/mdrun.out files. What I see in the regresiontests output is
that the forces/energies results are simply not correct; some tests simply
crash after 1-2 steps, but others do complete (like the nbnxn-free-energy/)
and the short-range energies a clearly far off.

I suggest to try to check if there may be hardware issue:

- run this memory testing tool:
git clone https://github.com/ComputationalRadiationPhysics/cuda_memtest.git
cd cuda_memtest
make cuda_memtest CFLAGS='-arch sm_30 -DSM_20 -O3 -DENABLE_NVML=0'
./cuda_memtest

- compile and run the gpu-burn tool:
git clone https://github.com/wilicc/gpu-burn
cd gpu-burn
make
then run
gpu-burn 300
to test for 5 minutes.

--
Szilárd


On Thu, Mar 28, 2019 at 3:46 PM Tafelmeier, Stefanie <
stefanie.tafelme...@zae-bayern.de> wrote:

> Hi Szilárd,
>
> Thanks again!
>
> Regarding the test:
>   -ntmpi 1 -ntomp 22 -pin on -pinstride 1:  2 out of 5 run
> https://it-service.zae-bayern.de/Team/index.php/s/XEQrYqq4pikGmMy  /
> https://it-service.zae-bayern.de/Team/index.php/s/YBdKKJ9c7zQpEg9
> Including:
>   -nsteps 0 -nb gpu -pme cpu -bonded cpu:   0 run
> https://it-service.zae-bayern.de/Team/index.php/s/YiByc7iXW5AW9ZX
>   -nsteps 0 -nb gpu -pme gpu -bonded cpu:   2 out of 5 run
> https://it-service.zae-bayern.de/Team/index.php/s/JNPXQnEgYtTAxGj   /
> https://it-service.zae-bayern.de/Team/index.php/s/6aq6BQwwbBELqWe
>   -nsteps 0 -nb gpu -pme gpu -bonded gpu:   0 run
> https://it-service.zae-bayern.de/Team/index.php/s/yj4RAqPMFsDNgTc
>
> Including:
>   -ntmpi 1 -ntomp 22 -pin on -pinstride 2:  1 out of 5 run
> https://it-service.zae-bayern.de/Team/index.php/s/q5jHbdJ2EygtDaQ  /
> https://it-service.zae-bayern.de/Team/index.php/s/sRPccwHRxojW9J8
>   -nsteps 0 -nb gpu -pme cpu -bonded cpu:   0 run
> https://it-service.zae-bayern.de/Team/index.php/s/GdKk5N68CY7BGxJ
>   -nsteps 0 -nb gpu -pme gpu -bonded cpu:   1 out of 5 run
> https://it-service.zae-bayern.de/Team/index.php/s/orwzKJMampWwDo5  /
> https://it-service.zae-bayern.de/Team/index.php/s/JXApT4tFtxQWxG6
>   -nsteps 0 -nb gpu -pme gpu -bonded gpu:   0 run
> https://it-service.zae-bayern.de/Team/index.php/s/8YKK7Zxax22RfGQ
>
> Including:
>   -ntmpi 1 -ntomp 22 -pin on -pinstride 4:  1 out of 5 run
> https://it-service.zae-bayern.de/Team/index.php/s/szZjzaxmwfimrgB  /
> https://it-service.zae-bayern.de/Team/index.php/s/QdTd2an9dbE9BSt
>   -nsteps 0 -nb gpu -pme cpu -bonded cpu:   3 out of 5 run
> https://it-service.zae-bayern.de/Team/index.php/s/DPoqKrgcWfF5PKM  /
> https://it-service.zae-bayern.de/Team/index.php/s/3NbsGHtCPsf7zFS
>   -nsteps 0 -nb gpu -pme gpu -bonded cpu:   3 out of 5 run
> https://it-service.zae-bayern.de/Team/index.php/s/WqP4tXjrR8i3455  /
> https://it-service.zae-bayern.de/Team/index.php/s/DACGc86xxKR6pWs
>   -nsteps 0 -nb gpu -pme gpu -bonded gpu:   0 run
> https://it-service.zae-bayern.de/Team/index.php/s/3nKdwA28KySLEdB
>
>
> Regarding the regressiontest:
> Here is the link to the tarball:
> https://it-service.zae-bayern.de/Team/index.php/s/mMyt3MPEfRrn8Ge
>
>
> Thanks again for all your support and fingers crossed!
>
> Best wishes,
> Steffi
>
>
>
>
>
> -Ursprüngliche Nachricht-
> Von: gromacs.org_gmx-users-boun...@maillist.sys.kth.se [mailto:
> gromacs.org_gmx-users-boun...@maillist.sys.kth.se] Im Auftrag von Szilárd
> Páll
> Gesendet: Mittwoch, 27. März 2019 20:27
> An: Discussion list for GROMACS users
> Betreff: Re: [gmx-users] WG: WG: Issue with CUDA and gromacs
>
> Hi Steffi,
>
> On Wed, Mar 27, 2019 at 1:08 PM Tafelmeier, Stefanie <
> stefanie.tafelme...@zae-bayern.de> wrote:
>
> > Hi Szilárd,
> >
> > thanks again!
> > Here are the links for the log files, that didn't run:
> > Old patch:
> >  -ntmpi 1 -ntomp 22 -pin on -pinstride 1:none ran*
> > https://it-service.zae-bayern.de/Team/index.php/s/b4AYiMCoHeNgJH3
> >  -ntmpi 1 -ntomp 22 -pin on -pinstride 2:none ran*
> > https://it-service.zae-bayern.de/Team/index.php/s/JEP2iwFFZCebZLF
> >  -ntmpi 1 -ntomp 22 -pin on -pinstride 4:one out of 5 ran
> > https://it-service.zae-bayern.de/Team/index.php/s/apra2zS7FHdqDQy
> >
> > New patch:
> >  -ntmpi 1 -ntomp 22 -pin on -pinstride 1:none ran*
> > https://it-service.zae-bayern.de/Team/index.php/s/jAD52jBgNddrS3w
> >  -ntmpi 1 -ntomp 22 -pin on -pinstride 2:none ran*
> > https://it-service.zae-bayern.de/Team/index.php/s/bcRjtz7r9NekzKB
> >  -ntmpi 1 -ntomp 22 -pin on -pinstride 4:none ran*
> > https://it-service.zae-bayern.de/Team/index.php/s/b3zp8DNztjE6ssF
> >
>
> This still doesn't tell much more unfortunately.
>
> Two more things to try (can be combined)
> - please set build with setting first
> cmake . -DCMAKE_BUILD_TYPE=RelWithAssert
> this may give us some extra debugging information during runs
> - please use this patch now -- it will print some 

[gmx-users] PME constant shifts.

[Sergio Perez]

Dear gmx comunity,

I am trying to calculate the electrostatic interaction of my system for
force field development. My system consists of 7 positively charged
particles interacting (system A) with another positively charged particle
(system B). I calculate the electrostatic interaction by doing reruns of
the trajectory in which some of them have been removed.

E_int(electrostatic) =
E_AB(coul-SR)+E_AB(coul-LR)-E_A(coul-SR)-E_A(coul-LR)-E_B(coul-SR)-E_B(coul-LR)

And to my surprise I get negative energies. I know the PME method has its
energy shifted:

"With the Verlet cut-off scheme, the PME direct space potential is shifted
by a constant such that
the potential is zero at the cut-off. This shift is small and since the net
system charge is close to
zero, the total shift is very small, unlike in the case of the
Lennard-Jones potential where all shifts
add up. We apply the shift anyhow, such that the potential is the exact
integral of the force."

Since for the system with just one charge I get a possitive E_1q(coul-LR),
I imagine the LR-coul term is shifted (eventhough I have to use the group
cut-off scheme and not the verlet). The problem is why don't all these
shifts cancel? Is there a way to not add the shifts?

I leave bellow the .mdp parameters. Note that for the VDW part I need to
use user-tables and therefore I need to use the group scheme.

Thank you very much!
Sergio Pérez-Conesa


; Integrator stuff
integrator = md-vv   ; steep or md
dt = 0.001
nsteps = 3 ; -1 no maximum
init-step = 30305000 ;
ld-seed = -2018962629
continuation = yes
emtol   = 10.0
emstep  = 0.01

; neighbour search
cutoff-scheme = group
nst-list = 10
verlet-buffer-tolerance = 0.005
ns-type = grid
; Pbc
pbc = xyz ; xyz or xy

; Vdw
rvdw = 1.4
rlist = 1.4
vdwtype = user  ; PME or User to look for table /user
DispCorr = No   ; corrections to energy and pressure or No

;Electrostatic
coulombtype = PME   ; User if look for table
rcoulomb = 1.4
pme-order =  4
fourierspacing = 0.2
ewald-geometry = 3d ; 3d or 3dc

;IMD-group=

; T and P
tcoupl =no  ; none, nose-hoover, v-rescale
nh-chain-length = 10
nsttcouple = 10 ; 10 if not used
pcoupl = no ;
pcoupltype = semiisotropic ;
compressibility = 4.5e-5 4.5e-5  ; isothermal compressibility of
water, bar^-1
refcoord-scaling = com
tau-p = 5.0  ; ps
ref-p = 1.0 1.0 ; bar

; Constraints
constraints = none  ; constraints  distances
constraint-algorithm = LINCS ; or SHAKE
lincs_iter  = 1
lincs_order = 4

;Generate velocities
gen_vel = no; assign velocities from Maxwell
distribution
gen_temp= 300   ; temperature for Maxwell distribution
gen_seed= 180611; generate a random seed

; OUTPUT
; Output control
nstenergy   = 1 ; save energies
nstlog  = 10 ; update log file every
;nstxout-compressed  = 500  ; save compressed coordinates every 10.0 ps
;nstxout = 20
;nstvout = 20
;nstfout = 20

; COM removal . There are more options like removing from groups or the
frequency
;comm-grps = System
comm-mode = none; linear, angular (removes linear and angular motion)
or none
nstcomm = 0  ; frequency of removal of com motion
freezegrps =
energygrp-excl =
;freezedim = Y Y Y
nstcalcenergy =
-- 
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] Energy Conservation at the Beginning of a Production Run

[Mark Abraham]

Hi,

Off topic for this discussion, but what would people think if we changed
things around so that grompp would do that constraining, before writing the
tpr file?

This will mean that the configuration in the tpr does not always match the
configuration that was given to grompp, but instead matches the rest of the
simulation. Thus such an energy jump would not occur. And it may mean that
broken structures are found at grompp time rather than later at step 0 when
doing mdrun on the cluster.

Mark

On Fri., 29 Mar. 2019, 02:17 Justin Lemkul,  wrote:

>
>
> On 3/28/19 6:02 PM, Kruse, Luke E.(MU-Student) wrote:
> > In the production run I have constraints on h-bonds and I do not have
> any constraints in the minimization.
>
> This would account for a sudden change in energy. At step 0, your H
> positions will immediately change.
>
> -Justin
>
> > 
> > From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <
> gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of David van
> der Spoel 
> > Sent: Thursday, March 28, 2019 4:50:38 PM
> > To: gmx-us...@gromacs.org
> > Subject: Re: [gmx-users] Energy Conservation at the Beginning of a
> Production Run
> >
> > Den 2019-03-28 kl. 21:19, skrev Kruse, Luke E.(MU-Student):
> >> In the mdout.mdp file corresponding to this run I have the line
> >>
> >> ; Do not constrain the start configuration
> >> continuation = no
> >>
> >>
> >> Do I need to change this to yes?
> > Maybe, try a 100 step run where you save energies at each time step. Did
> > you change anything else between minimization and production, e.g.
> > constraint settings?
> >>
> >> The plot of energy vs time looks like an under damped response to a set
> point change in control theory - slight overshooting after the first few
> steps, then settling down to an average value nicely. This is what is
> leading me to believe that there is a large, non conservative force near
> the beginning that becomes smaller, and conservative.
> >>
> >> Time (ps)
> >>
> >>   0.00  -1851396.981075
> >>   0.10  -1851388.994856
> >>   0.20  -1851390.001591
> >>   0.30  -1851389.436699
> >>   0.40  -1851389.760231
> >>   0.50  -1851389.726310
> >>   0.60  -1851389.808428
> >>   0.70  -1851389.907829
> >>   0.80  -1851389.840644
> >>   0.90  -1851389.967373
> >>   1.00  -1851390.061502
> >>   1.10  -1851390.057500
> >>   1.20  -1851390.086172
> >>   1.30  -1851390.042747
> >>   1.40  -1851390.060210
> >>   1.50  -1851389.998599
> >>   1.60  -1851390.120319
> >>   1.70  -1851390.141846
> >>   1.80  -1851390.131003
> >>   1.90  -1851390.170220
> >>   2.00  -1851390.141617
> >>   2.10  -1851390.135779
> >>   2.20  -1851390.255565
> >>   2.30  -1851390.153947
> >>   2.40  -1851390.238681
> >>   2.50  -1851390.144512
> >>   2.60  -1851390.255200
> >>   2.633000  -1851390.178012
> >>
> >>
> >> 
> >> From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <
> gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of David van
> der Spoel 
> >> Sent: Thursday, March 28, 2019 3:10:05 PM
> >> To: gmx-us...@gromacs.org
> >> Subject: Re: [gmx-users] Energy Conservation at the Beginning of a
> Production Run
> >>
> >> Den 2019-03-28 kl. 20:53, skrev Kruse, Luke E.(MU-Student):
> >>> Hello gromacs users,
> >>>
> >>>
> >>> I am trying to simulate a peptide amphiphile with the CHARMM27 force
> field. To do this I have had to specify an additional bond type and bond,
> angle, dihedral, etc. parameters in the .itp files of the force field.
> Then, to check if I had done this correctly, I minimized the system, and
> ran a production run to generate an NVE ensemble, so that I could make sure
> the system was conserving energy appropriately. After looking at the energy
> vs time plot (produces with the gmx energy command), however, the system
> jumps from an initial energy, up ~20 kJ per mole and then conserves energy
> for the most part (slight drifting but seems tolerable relative to the
> initial discontinuity). Is this discontinuity a normal happenstance or a
> result of bad minimization?
> >>>
> >> Is that at step 0 in the simulation?
> >>
> >> Please check options like
> >> ; Do not constrain the start configuration
> >> continuation = no
> >>
> >>
> >>> Thank you!
> >>>
> >>> Luke
> >>>
> >>
> >> --
> >> David van der Spoel, Ph.D., Professor of Biology
> >> Head of Department, Cell & Molecular Biology, Uppsala University.
> >> Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
> >> http://www.icm.uu.se
> >> --
> >> 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 

Re: [gmx-users] Coarse-grained Protein-ligand simulations

[Mac Kevin Braza]

Hello everyone,

I am simulating a coarse-grained model of a membrane protein (GPCR) in
lipid bilayer and an all-atom ligand octopamine. I build the protein,
solutes, and membrane in the web server CHARMM-GUI. While, I added the
ligand to the protein complex manually using the same coordinates of the
coarse-grained protein model.

I used the GROMACS input files from the output of CHARMM-GUI to simulate
the system. I include the LIGAND.ITP (from the PRODRG Server) to the
system.top and added the atom indexes in the index.ndx file.

However, when I proceed with the second part of equilibration, the
following errors occurred.

*Command line*:
  gmx grompp -f step6.2_equilibration.mdp -o step6.2_equilibration.tpr -c
step6.1_equilibration.gro -p system.top -n index.ndx

Setting the LD random seed to 1722366284
Generated 2391 of the 4656 non-bonded parameter combinations
Excluding 1 bonded neighbours molecule type 'PROA_P'
Excluding 1 bonded neighbours molecule type 'POPC'
Excluding 1 bonded neighbours molecule type 'W'
Excluding 1 bonded neighbours molecule type 'NA'
Excluding 1 bonded neighbours molecule type 'CL'
Excluding 3 bonded neighbours molecule type 'LIG'
Velocities were taken from a Maxwell distribution at 303.15 K
Removing all charge groups because cutoff-scheme=Verlet

---
Program gmx grompp, VERSION 5.1.4
Source code file: /home/gromacs-5.1.4/src/gromacs/gmxpreprocess/readir.c,
line: 2690

Fatal error:
20 atoms are not part of any of the T-Coupling groups
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---

The 20 atoms described the ligand I placed inside the protein-membrane
complex. I want to know if where can this error originate and how can we
fix them?

We will appreciate any help you will give. Thank you very much!

Best regards,
Mac Kevin E. Braza

On Fri, Mar 29, 2019 at 2:28 PM Mac Kevin Braza  wrote:

> I am simulating a coarse-grained model of a octopamine receptor in lipid
> bilayer and an all-atom ligand octopamine. I build the protein, solutes,
> and membrane in the web server CHARMM-GUI. While, I added the ligand to the
> protein complex manually using the same coordinates of the coarse-grained
> protein model.
>
> I used the GROMACS input files from the output of CHARMM-GUI to simulate
> the system. I include the LIGAND.ITP (from the PRODRG Server) to the
> system.top and added the atom indexes in the index.ndx file. I attached
> here the two files.
>
> However, when I proceed with the second part of equilibration, the
> following errors occurred.
>
> *Command line*:
>   gmx grompp -f step6.2_equilibration.mdp -o step6.2_equilibration.tpr -c
> step6.1_equilibration.gro -p system.top -n index.ndx
>
> Setting the LD random seed to 1722366284
> Generated 2391 of the 4656 non-bonded parameter combinations
> Excluding 1 bonded neighbours molecule type 'PROA_P'
> Excluding 1 bonded neighbours molecule type 'POPC'
> Excluding 1 bonded neighbours molecule type 'W'
> Excluding 1 bonded neighbours molecule type 'NA'
> Excluding 1 bonded neighbours molecule type 'CL'
> Excluding 3 bonded neighbours molecule type 'LIG'
> Velocities were taken from a Maxwell distribution at 303.15 K
> Removing all charge groups because cutoff-scheme=Verlet
>
> ---
> Program gmx grompp, VERSION 5.1.4
> Source code file:
> /home/minco/gromacs-5.1.4/src/gromacs/gmxpreprocess/readir.c, line: 2690
>
> Fatal error:
> 20 atoms are not part of any of the T-Coupling groups
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
> ---
>
> The 20 atoms described the ligand I placed inside the protein-membrane
> complex. I want to know if where can this error originate and how can we
> fix them?
>
> We will appreciate any help you will give. Thank you very much!
>
> Best regards,
> Mac Kevin E. Braza
>
-- 
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.