normal Clarifiaction about chlrophyll a molecule Coarsed grained topology file

  • estebanp
  • estebanp's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
7 years 7 months ago #3722 by estebanp
Dear MARTINI users,
I have a problem to include a positions restraints file in the topology file (gromacs). My system is 8x8x1 DPPC lipids coarse grained (I took 1 single molecule from the dppc_water.gro, erasing the rest and the water) and construct the layer with genconf.
I center the system (0,0,0) in a bigger box with sizes 10.5532 11.68792 6.105482
I make genrestr for generate the posre.itp, selecting the only group thet exists, DPPC. Then I make my topology like:

&&&&&&&&&&&&

#include "/home/pedrueza/monolayer/dppc_single/cg/martini_v2.0.itp"
#include "/home/pedrueza/monolayer/dppc_single/cg/martini_v2.0_lipids.itp"

; Include Position restraint file
#ifdef POSRES
#include "/home/pedrueza/monolayer/dppc_single/cg/posre.itp"
#endif

[ system ]
DPPC MONOLAYER

[ molecules ]
DPPC 64

&&&&&&&&&

Then I minimaze the energy. First, without “define = -DPOSRES ; position restrain the lipids” in the .mdp file, and everything is OK (surprisingly, in VMD the molecules seems that they don`t move!!!). But when I activate the -DPOSRES and run grompp, I get the following error:


Fatal error:
[ file posre.itp, line 29 ]:
Atom index (25) in position_restraints out of bounds (1-24).
This probably means that you have inserted topology section "position_restraints"
in a part belonging to a different molecule than you intended to.
In that case move the "position_restraints" section to the right molecule.

One solution I probe is to insert the “; Include Position restraint file
“ lines into the martini_v2.0_lipids.itp, at the end of DPPC definition, by this way:


&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

;;;;;; DIPALMITOYL PHOSPHATIDYLCHOLINE
;
; in general models PCs with saturated tail lengths C15-18

[moleculetype]
; molname nrexcl
DPPC 1

[atoms]
; id type resnr residu atom cgnr charge
1 Q0 1 DPPC NC3 1 1.0
2 Qa 1 DPPC PO4 2 -1.0
3 Na 1 DPPC GL1 3 0
4 Na 1 DPPC GL2 4 0
5 C1 1 DPPC C1A 5 0
6 C1 1 DPPC C2A 6 0
7 C1 1 DPPC C3A 7 0
8 C1 1 DPPC C4A 8 0
9 C1 1 DPPC C1B 9 0
10 C1 1 DPPC C2B 10 0
11 C1 1 DPPC C3B 11 0
12 C1 1 DPPC C4B 12 0

[bonds]
; i j funct length force.c.
1 2 1 0.470 1250
2 3 1 0.470 1250
3 4 1 0.370 1250
3 5 1 0.470 1250
5 6 1 0.470 1250
6 7 1 0.470 1250
7 8 1 0.470 1250
4 9 1 0.470 1250
9 10 1 0.470 1250
10 11 1 0.470 1250
11 12 1 0.470 1250

[angles]
; i j k funct angle force.c.
2 3 4 2 120.0 25.0
2 3 5 2 180.0 25.0
3 5 6 2 180.0 25.0
5 6 7 2 180.0 25.0
6 7 8 2 180.0 25.0
4 9 10 2 180.0 25.0
9 10 11 2 180.0 25.0
10 11 12 2 180.0 25.0

; Include Position restraint file
#ifdef POSRES
#include "/home/pedrueza/monolayer/dppc_single/cg/posre.itp"
#endif

;;;;;; DIHEXANOYL PHOSPHATIDYLCHOLINE
;
; in general models PCs with saturated tail lengths C8-11

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


but again I obtain the same error, this time with different numbers:

Fatal error:
[ file posre.itp, line 17 ]:
Atom index (13) in position_restraints out of bounds (1-12).
This probably means that you have inserted topology section "position_restraints"
in a part belonging to a different molecule than you intended to.
In that case move the "position_restraints" section to the right molecule.


I understand that the numeration must be incorrect, but it seems correct when I see the files. The make_ndx command generates the same type of file that following the genrestr options. Maybe is a problem of isolated one DPPC molecule?

Can anyone help me?. I will be gratefull with any help that you can provide me.

P.D.:
the header and end of my .gro file is:

DPPC BILAYER
768
1DPPC NC3 1 -4.349 -4.676 -0.841 -0.0814 0.0799 0.0504
1DPPC PO4 2 -4.690 -4.928 -0.916 -0.0933 -0.0920 0.2849
1DPPC GL1 3 -4.708 -5.098 -0.588 -0.2760 0.0809 0.1057
1DPPC GL2 4 -4.593 -5.437 -0.612 -0.0911 -0.0996 0.2134
1DPPC C1A 5 -4.743 -4.894 -0.210 -0.1749 -0.2140 0.2303
1DPPC C2A 6 -4.467 -4.942 0.138 -0.1136 -0.2962 -0.0951
1DPPC C3A 7 -4.334 -4.834 0.580 -0.0554 -0.0277 -0.0337
1DPPC C4A 8 -4.039 -4.484 0.722 0.0471 -0.2678 0.1634
1DPPC C1B 9 -4.706 -5.498 -0.176 0.0266 -0.1530 -0.1297
1DPPC C2B 10 -4.605 -5.433 0.309 -0.0409 0.1791 0.0088
1DPPC C3B 11 -4.907 -5.478 0.685 0.2918 -0.0020 -0.1119
1DPPC C4B 12 -5.258 -5.659 0.903 0.1992 -0.0507 0.0256
2DPPC NC3 13 -4.349 -3.215 -0.841 -0.0814 0.0799 0.0504
2DPPC PO4 14 -4.690 -3.467 -0.916 -0.0933 -0.0920 0.2849

(...)


63DPPC C2B 754 4.629 3.333 0.309 -0.0409 0.1791 0.0088
63DPPC C3B 755 4.327 3.288 0.685 0.2918 -0.0020 -0.1119
63DPPC C4B 756 3.976 3.107 0.903 0.1992 -0.0507 0.0256
64DPPC NC3 757 4.885 5.551 -0.841 -0.0814 0.0799 0.0504
64DPPC PO4 758 4.544 5.299 -0.916 -0.0933 -0.0920 0.2849
64DPPC GL1 759 4.526 5.129 -0.588 -0.2760 0.0809 0.1057
64DPPC GL2 760 4.641 4.790 -0.612 -0.0911 -0.0996 0.2134
64DPPC C1A 761 4.491 5.333 -0.210 -0.1749 -0.2140 0.2303
64DPPC C2A 762 4.767 5.285 0.138 -0.1136 -0.2962 -0.0951
64DPPC C3A 763 4.900 5.393 0.580 -0.0554 -0.0277 -0.0337
64DPPC C4A 764 5.195 5.743 0.722 0.0471 -0.2678 0.1634
64DPPC C1B 765 4.528 4.729 -0.176 0.0266 -0.1530 -0.1297
64DPPC C2B 766 4.629 4.794 0.309 -0.0409 0.1791 0.0088
64DPPC C3B 767 4.327 4.749 0.685 0.2918 -0.0020 -0.1119
64DPPC C4B 768 3.976 4.568 0.903 0.1992 -0.0507 0.0256
10.55320 11.68792 6.10548


and the posre.itp generated with genrestr

; position restraints for DPPC of DPPC BILAYER

[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
2 1 1000 1000 1000
3 1 1000 1000 1000
4 1 1000 1000 1000
5 1 1000 1000 1000
6 1 1000 1000 1000
7 1 1000 1000 1000

(...)

761 1 1000 1000 1000
762 1 1000 1000 1000
763 1 1000 1000 1000
764 1 1000 1000 1000
765 1 1000 1000 1000
766 1 1000 1000 1000
767 1 1000 1000 1000
768 1 1000 1000 1000

Please Log in or Create an account to join the conversation.

More
7 years 7 months ago #3725 by jaakko
Hi,

the position restraint file can only have force constants for each atom/bead in the itp. So if you include it in the molecule itp (second version you tried) it should only have indices up to 12. Remove the rest of the lines from your posre.itp and you should be fine.

Btw. the reason why the first thing you tried complains about index 25 is that the last lipid defined in martini_v2.0_lipids.itp has 24 beads. Adding the #ifdef POSRES after the #includes in the .top file is basically adding it to the last molecule definition in the last #include. You should think all the #includes as copy pasting the whole file to that place in the .top file.

Please Log in or Create an account to join the conversation.

  • estebanp
  • estebanp's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
7 years 7 months ago #3726 by estebanp
Thank you very much, jaako. I sucessfully solve the problem deleting the other lines, as you said. :)
Cheers,
Esteban

Please Log in or Create an account to join the conversation.

More
9 months 3 weeks ago #8829 by saini
I am also getting same error. I want to apply position restraint on lipid.
; molname nrexcl
PPG 1

[ atoms ]
; id type resnr residu atom cgnr charge
1 P4 1 PPG GL0 1 0 72.0000
2 Qa 1 PPG PO4 2 -1 72.0000
3 Na 1 PPG GL1 3 0 72.0000
4 Na 1 PPG GL2 4 0 72.0000
5 C1 1 PPG C1A 5 0 72.0000
6 C1 1 PPG C2A 6 0 72.0000
7 C1 1 PPG C3A 7 0 72.0000
8 C1 1 PPG C4A 8 0 72.0000
9 C3 1 PPG D1B 9 0 72.0000
10 C1 1 PPG C2B 10 0 72.0000
11 C1 1 PPG C3B 11 0 72.0000
12 C1 1 PPG C4B 12 0 72.0000

[ bonds ]
; i j funct length force.c.
1 2 1 0.470 1250
2 3 1 0.470 1250
3 4 1 0.370 1250
3 5 1 0.470 1250
5 6 1 0.470 1250
6 7 1 0.470 1250
7 8 1 0.470 1250
4 9 1 0.470 1250
9 10 1 0.470 1250
10 11 1 0.470 1250
11 12 1 0.470 1250
angles ]
; i j k funct angle force.c.
2 3 4 2 120.0 25.0
2 3 5 2 180.0 25.0
3 5 6 2 180.0 25.0
5 6 7 2 180.0 25.0
6 7 8 2 180.0 25.0
4 9 10 2 180.0 45.0
9 10 11 2 180.0 25.0
10 11 12 2 180.0 25.0

#ifdef POSRES_LIPID
#include "ppg_cg_pr.itp"
#endif


;phosphatidyl-glycerole (PG), with a trans-3-hexadecenoic acid tail (sn1) and an α-linolenoyl tail (sn2)

[ moleculetype ]
; molname nrexcl
PPT 1

[ atoms ]
; id type resnr residu atom cgnr charge
1 P4 1 PPT GL0 1 0 72.0000

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

I have generate position restraint on by using gmx genrestr command

it looks like this
; position restraints for PPG of Protein DPPC BILAYER

[ position_restraints ]
; i funct fcx fcy fcz
2377 1 1000 1000 1000
2378 1 1000 1000 1000
2379 1 1000 1000 1000
2380 1 1000 1000 1000
2381 1 1000 1000 1000
2382 1 1000 1000 1000
2383 1 1000 1000 1000
2384 1 1000 1000 1000
2385 1 1000 1000 1000
2386 1 1000 1000 1000
2387 1 1000 1000 1000
2388 1 1000 1000 1000
2389 1 1000 1000 1000
2390 1 1000 1000 1000
2391 1 1000 1000 1000
2392 1 1000 1000 1000
2393 1 1000 1000 1000
2394 1 1000 1000 1000
2395 1 1000 1000 1000

........................
17185 1 1000 1000 1000
17186 1 1000 1000 1000
17187 1 1000 1000 1000
17188 1 1000 1000 1000
17189 1 1000 1000 1000
17190 1 1000 1000 1000
17191 1 1000 1000 1000
17192 1 1000 1000 1000
17193 1 1000 1000 1000
17194 1 1000 1000 1000
17195 1 1000 1000 1000
after running grompp command for minimization i got the error
ERROR 1 : Atom index (2377) in position_restraints out of bounds (1-12). This probably means that you have inserted topology section "position_restraints" in a part belonging to a different molecule than you intended to. In that case move the "position_restraints" section to the right molecule Please suggest me, how can i solve it?[file ppg_cg_pr.itp, line 5]:
Atom index (2377) in position_restraints out of bounds (1-12).
This probably means that you have inserted topology section
"position_restraints"
in a part belonging to a different molecule than you intended to.
In that case move the "position_restraints" section to the right molecule

Please suggest me, how can i solve it?

Please Log in or Create an account to join the conversation.

More
9 months 3 weeks ago #8830 by vainikka
From a first glance it looks you made the same mistake as the original poster. The postiion restraints are supposed to be for each bead in the itp file, not in the gro file. Since you do not have 2377 beads defined in the itp, GROMACS gives you an error.

Please Log in or Create an account to join the conversation.

More
9 months 3 weeks ago #8831 by saini
Thank you so much for kindly reply. I got it and successfully solve the problem by deleting other lines.
now my posre.itp looks like this.

; position restraints for PPG of Protein DPPC BILAYER

[ position_restraints ]
; i funct fcx fcy fcz
2 1 1000 1000 1000

Please Log in or Create an account to join the conversation.

More
9 months 3 weeks ago #8832 by saini
Dear MARTINI users,
I am trying to insert protein_pigment complex by using g_membed command. g_membed.mdp file is

title = Martini
cpp = /usr/bin/cpp
define = -DPOSRES -DPOSRES_LIPID
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 10000
init_step = 0
comm-mode = Linear
nstcomm = 100
comm-grps =

; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 1000
emstep = 0.01
; Max number of iterations in relax_shells
niter = 20
; Step size (ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 1000
nbfgscorr = 10

; TEST PARTICLE INSERTION OPTIONS
rtpi = 0.05

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 5000
nstvout = 5000
nstfout = 5000
; Output frequency for energies to log file and energy file
nstlog = 5000
nstenergy = 100
; Output frequency and precision for xtc file
nstxtcout = 1000
xtc-precision = 100
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps =Protein_CLA_CHL_LUT

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 10
; ns algorithm (simple or grid)
ns-type = Grid
; Periodic boundary conditions: xyz, no, xy
pbc = xyz
periodic_molecules = no
; nblist cut-off
rlist = 1.2

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = Cut-off
rcoulomb-switch = 0
rcoulomb = 1.2
; Relative dielectric constant for the medium and the reaction field
epsilon_r = 1
epsilon_rf = 1
; Method for doing Van der Waals
vdw-type = Cut-off
; cut-off lengths
rvdw-switch = 0
rvdw = 1.2
cutoff-scheme = group
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = No
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; Seperate tables between energy group pairs
energygrp_table =
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = v-rescale
; Groups to couple separately
tc-grps = Protein_lipids sol
; Time constant (ps) and reference temperature (K)
tau-t = 0.1 0.1
ref-t = 293 293
; Pressure coupling
Pcoupl = no
Pcoupltype = semiisotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1
compressibility = 5.3e-05 5.3e-05
ref-p = 1.01325 1.01325
; Scaling of reference coordinates, No, All or COM
refcoord_scaling = No
; Random seed for Andersen thermostat
andersen_seed = 815131

; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel = yes
gen-temp = 323
gen-seed = 173529

; OPTIONS FOR BONDS
constraints = none ;all-bonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
continuation = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 0.0001
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 30
; Convert harmonic bonds to morse potentials
morse = no
; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl = Protein_CLA_CHL_LUT Protein_CLA_CHL_LUT


; Non-equilibrium MD stuff
acc-grps =
accelerate =
freezegrps =Protein_CLA_CHL_LUT
freezedim =Y Y Y



and membed.dat file is

nxy = 1000 ;Number of iteration for the xy dimension
nz = 0 ; Number of iterations for the z dimension
xyinit = 0.50000 ; Resize factor for the protein in the xy dimension before starting embedding
xyend = 1.000000 ; Final resize factor in the xy dimension
zinit = 1.000000 ;Resize factor for the protein in the z dimension before starting embedding
zend = 1.000000 ;Final resize faction in the z dimension
maxwarn = 0
~
but after gmx mdrun i am getting LINCS warning.
Please suggest me what changes should i have to do in .mdp and .dat file.

Please Log in or Create an account to join the conversation.

More
9 months 3 weeks ago #8835 by vainikka
Is this a CG system? If yes, then I would not necessarily use the GWDG example mdp file since it is meant for atomistic systems, at least as far as I can tell.

If you are working with a Martini membrane and protein, then I would recommend using insane to build such systems.

Please Log in or Create an account to join the conversation.

More
9 months 3 weeks ago #8837 by saini
First of all, thank you for kindly reply.
Yes, it is a CG system. I build a protein membrane system by using insane.py script.
by using following command
python insane.py -f lhcii_em.gro -o system.gro -p system.top -pbc cubic -box 25,25,9 -l FPGG:14 -l DFGG:72 -l FPMG:14 -l DFMG:100 -l JFPG:14 -l FPSG:42 -l JPPG:14 -u FPGG:14 -u DFGG:72 -u FPMG:14 -u DFMG:100 -u JFPG:14 -u FPSG:42 -u JPPG:14 -salt 0.1 -center -sol W

Then i did enenrgy minimization after that i am doing nvt equilibriation by using the command
Fatal error:
An atom moved too far between two domain decomposition steps
This usually means that your system is not well equilibrated

i am using this mdp file is
; VARIOUS PREPROCESSING OPTIONS =
title = Martini
;define = -DPOSRES
; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit = 0.0
dt = 0.002
nsteps = 100000 ;= 200ps
; number of steps for center of mass motion removal =
nstcomm = 100
comm-grps = Protein_lipids sol

; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 0
nstvout = 0
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 5000
nstenergy = 5000
; Output frequency and precision for xtc file =
nstxtcout = 5000
xtc_precision = 5000
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
xtc-grps = Protein_lipids sol
; Selection of energy groups =
energygrps = Protein_lipids sol

; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist = 10
; ns algorithm (simple or grid) =
ns_type = grid
; Periodic boundary conditions: xyz or none =
pbc = xyz
; nblist cut-off =
rlist = 1.2

; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = Cut-off
rcoulomb_switch = 0.0
rcoulomb = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r = 15
; Method for doing Van der Waals =
vdw_type = Cut-off
; cut-off lengths =
rvdw_switch = 0.9
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = v-rescale
; Groups to couple separately =
tc-grps = Protein_lipids sol
; Time constant (ps) and reference temperature (K) =
tau_t = 2.0 2.0
ref_t = 293 293

; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel = yes
gen_temp = 105
gen_seed = 548628
; OPTIONS FOR BONDS =
constraints = none
; Type of constraint algorithm =
constraint_algorithm = Lincs
; Highest order in the expansion of the constraint coupling matrix =
lincs_order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs_warnangle = 30



Please let me know, where i am doing mistake?

Please Log in or Create an account to join the conversation.

More
9 months 2 weeks ago #8838 by vainikka
Are you running the NVT simulation on all available cores? If yes, try running a few hundred ps on just a single core.

Please Log in or Create an account to join the conversation.

More
9 months 2 weeks ago #8839 by saini
yes, i am running NVT simulation after enenrgy minimization

I am using these command
gmx grompp -f martini_nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -v -deffnm nvt

then i tried by using -nt 1 option or -nt 8 option
gmx mdrun -v -deffnm nvt -nt 1

but i am getting error after 20 ps, after that i got LINCS warning.

Step 3059, time 6.118 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms -nan, max 1851868416.000000 (between atoms 872 and 874)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
660 661 90.1 0.5397 0.9609 0.2400
660 662 90.1 0.5089 0.8379 0.2400
661 662 90.2 0.2388 0.3047 0.2400
667 668 90.0 1.7553 3.0411 0.2000
690 691 43.4 0.2388 0.2344 0.2400
690 692 82.6 0.2482 0.2487 0.2400
691 692 38.6 0.2361 0.2767 0.2400
693 694 58.9 16.4492 40.7888 0.1950
872 873 104.5 0.2400 13.3778 0.2400
872 874 134.2 0.2400 444448416.0000 0.2400
873 874 76.1 0.2400 444448416.0000 0.2400
929 930 90.0 0.2200 4.9042 0.2200
930 931 90.0 0.2500 5.5867 0.2500
931 932 57.5 0.2800 0.2801 0.2800
932 929 90.0 0.2550 2.6561 0.2550
1919 1923 90.0 2.2908 3.1582 0.8000
2234 2238 65.0 0.8000 320.1271 0.8000
Wrote pdb files with previous and current coordinates


i have set costraint = none but i don't know why i am getting LINCS warning?

Please let me know, if you have any idea about it.

Please Log in or Create an account to join the conversation.

More
9 months 2 weeks ago #8840 by saini
In CG system
I am trying to run NVT simulation after energy minimization. it ran 10 ns sucessufully but after that i got the error
Fatal error:
Step 109300: The total potential energy is -nan, which is not finite. The LJ
and electrostatic contributions to the energy are 0 and -22524.4,
respectively. A non-finite potential energy can be caused by overlapping
interactions in bonded interactions or very large or Nan coordinate values.
Usually this is caused by a badly- or non-equilibrated initial configuration,
incorrect interactions or parameters in the topology.

what is the meaning of -nan?
Please suggest me something, how can i solve this error?

Please Log in or Create an account to join the conversation.

More
9 months 2 weeks ago #8841 by vainikka
NaN = Not a Number. Effectively your system blew up, and one or more terms has such a large value that GROMACS does not want to keep track of it anymore. You can check the GROMACS mailing list for more information.

The output you have copy-pasted here most likely contains the reason for your error: "Usually this is caused by a badly- or non-equilibrated initial configuration,
incorrect interactions or parameters in the topology".

Also see: manual.gromacs.org/documentation/2018/us...tml#system-diagnosis

Please Log in or Create an account to join the conversation.

More
9 months 2 weeks ago #8842 by saini
Ok, thank you so much for your reply but i did not got any solution on mailing list.

I am runnung the simulation further by using -cpt npt.cpt -append option but it stop every 1ns or 2ns run and gave same erreor. so i am using -cpi npt.cpt -append option everytime.

the command that i am using
gmx mdrun -v -deffnm npt -ntmpi 1 -ntomp 1 -pin on -tunepme no -nstlist 400 -nb gpu -pme gpu -cpi npt.cpt -append

it is going wrong, isn't it?

Please Log in or Create an account to join the conversation.

More
9 months 2 weeks ago #8849 by vainikka
Are you still using the same mdp file you posted here before? If yes, then what you are seeing is likely an error caused by faulty topology. Which system is it that you are simulating, and where did you get the itp file from?

Please Log in or Create an account to join the conversation.

More
9 months 2 weeks ago #8850 by saini
Thank you for your kindly conversation.
I am new Martini users so i am not understanding where i am doing mistake.

now i am doing NPT simulation after 500ps NVT run. but i am getting same error "potential energy -nan". my mdp file is

title = Martini
integrator = md
tinit = 0.0
dt = 0.004
nsteps = 250000000 ;1000ns
nstcomm = 100
comm-grps = cla_lipids sol
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 50000
nstenergy = 50000
nstxtcout = 50000
xtc_precision = 50000
xtc-grps = cla_lipids sol
cutoff-scheme = verlet
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.2
coulombtype = PME
pme-order = 4
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15
vdw_type = Cut-off
rvdw_switch = 0.9
rvdw = 1.2
DispCorr = No
tcoupl = Berendsen
tc-grps = cla_lipids sol
tau_t = 2.0 2.0
ref_t = 293 293
Pcoupl = Berendsen
Pcoupltype = semiisotropic
tau_p = 1
compressibility = 3e-4 3e-4
ref_p = 1.0 1.0
gen_vel = no
gen_temp = 105
gen_seed = 548628
continuation = yes
constraints = none
constraint_algorithm = Lincs
lincs_order = 4
lincs_warnangle = 30

My system is plant thylakoid bilayer and chlorophyll a (CLA) system.

I have equilibriated bilayer without CLA for 1 microsec NPT succesfully.
after that i added CLA molecules in equilibriated bilayer by "gmx editconf' command. Then did energy minimization sucussefully. but during equilibriation i got this error (potential energy ) after 5 or 10 ns run length.

I did not found CLA gro file at martini website. i have taken it from my freind.
the gro file
CLA
24
1CLA P5N 1 1.902 12.059 4.639 -0.2338 0.0458 0.0760
1CLA P1N_1 2 1.998 12.228 4.728 -0.2383 -0.2153 -0.0255
1CLA P1N_2 3 2.033 12.093 4.481 0.4977 -0.0777 0.1318
1CLA P1N_3 4 1.958 11.859 4.582 0.1736 0.0644 -0.4561
1CLA P1N_4 5 1.975 11.979 4.822 -0.1787 0.2577 -0.4076
1CLA SC3_1 6 1.999 12.184 4.962 -0.1448 -0.2055 0.3604
1CLA SC3_2 7 2.024 11.755 4.797 -0.3473 0.2860 0.0958
1CLA SC3_8 8 1.983 11.625 4.442 0.0938 -0.2695 0.0421
1CLA C2_1 9 2.085 11.493 4.684 -0.4512 -0.4162 0.1746
1CLA SC3_3 10 1.962 11.898 4.347 0.2372 -0.2194 -0.3866
1CLA SC3_4 11 2.184 12.196 4.288 0.0744 0.1325 -0.1576
1CLA C3_2 12 2.178 12.007 4.120 -0.2895 0.0432 0.3014
1CLA SC3_5 13 2.169 12.274 4.567 0.0133 0.1509 -0.0053
1CLA SC3_6 14 2.094 12.452 4.758 -0.4154 -0.0706 -0.3372
1CLA SC3_7 15 2.032 11.845 5.065 0.1887 0.3553 0.2883
1CLA SC5_1 16 2.088 12.032 5.208 -0.3525 0.1092 -0.3116
1CLA Na_1 17 1.940 12.303 5.176 0.1154 -0.4853 0.1600
1CLA Na_2 18 2.027 12.307 5.367 -0.0345 0.3314 -0.5669
1CLA C2_2 19 1.882 12.494 4.873 0.0969 0.0074 -0.1357
1CLA Na_3 20 1.749 12.670 4.871 0.0834 0.4260 -0.1194
1CLA C3_1 21 1.726 12.939 4.653 0.2145 0.0221 -0.0616
1CLA C1_1 22 1.916 13.084 4.219 0.4054 -0.0219 0.2319
1CLA C1_2 23 1.618 13.394 3.956 -0.0895 -0.0303 0.0300
1CLA C1_3 24 1.375 13.785 3.843 0.2168 0.1923 0.0360
10 10 10

and itp that i am using

[ moleculetype ]
CLA 3

[atoms]
;nr type resnr residue atom cgnr charge mass
1 P5N 1 CLA P5N 1 1.00000 24.305000
2 P1N 1 CLA P1N_1 2 -.25000 14.006700
3 P1N 1 CLA P1N_2 3 -.25000 14.006700
4 P1N 1 CLA P1N_3 4 -.25000 14.006700
5 P1N 1 CLA P1N_4 5 -.25000 14.006700
6 SC3 1 CLA SC3_1 6 0.000000 36.032999
7 SC3 1 CLA SC3_2 7 0.000000 37.040999
8 SC3 1 CLA SC3_8 8 0.000000 36.032999
9 C2 1 CLA C2_1 9 0.000000 27.046000
10 SC3 1 CLA SC3_3 10 0.000000 37.040999
11 SC3 1 CLA SC3_4 11 0.000000 39.056999
12 C3 1 CLA C3_2 12 0.000000 27.046000
13 SC3 1 CLA SC3_5 13 0.000000 37.040999
14 SC3 1 CLA SC3_6 14 0.000000 39.056999
15 SC3 1 CLA SC3_7 15 0.000000 39.056999
16 SC5 1 CLA SC5_1 16 0.000000 30.026401
17 Na 1 CLA Na_1 17 0.000000 41.029400
18 Na 1 CLA Na_2 18 0.000000 31.034400
19 C2 1 CLA C2_2 19 0.000000 28.054001
20 Na 1 CLA Na_3 20 0.000000 44.009800
21 C3 1 CLA C3_1 21 0.000000 69.127001
22 C1 1 CLA C1_1 22 0.000000 70.135001
23 C1 1 CLA C1_2 23 0.000000 70.135001
24 C1 1 CLA C1_3 24 0.000000 71.143001

[ bonds ]
1 2 1 0.214 61074.4 ; 1:bond1:1
1 3 1 0.212 61074.4 ; 1:bond2:1
1 4 1 0.212 54348.8 ; 1:bond3:1
1 5 1 0.21 59240.6 ; 1:bond4:1
2 6 1 0.237 91641.4 ; 1:bond5:1
2 13 1 0.246 70109.4 ; 1:bond6:1
3 13 1 0.244 85688.8 ; 1:bond7:1
3 10 1 0.242 90121.8 ; 1:bond8:1
4 10 1 0.24 119967.2 ; 1:bond9:1
4 7 1 0.243 124988.4 ; 1:bond10:1
5 7 1 0.241 121873.6 ; 1:bond11:1
5 6 1 0.25 121873.6 ; 1:bond12:1
6 15 1 0.3501 88505.6 ; 1:bond13:1
6 14 1 0.356 44661.4 ; 1:bond14:1
6 17 1 0.216 6780.88 ; 1:bond15:1
7 8 1 0.372 77226.2 ; 1:bond16:1
7 15 1 0.287 84689.6 ; 1:bond17:1
8 10 1 0.292 82460 ; 1:bond18:1
10 11 1 0.36 36158 ; 1:bond19:1
11 13 1 0.286 64931.6 ; 1:bond20:1
13 14 1 0.284 67572.8 ; 1:bond21:1
11 12 1 0.2738 23955.6 ; 1:bond22:1
8 9 1 0.292 44550.4 ; 1:bond23:1
15 16 1 0.261 32067.4 ; 1:bond24:1
16 17 1 0.308 11551.76 ; 1:bond25:1
17 18 1 0.207 83262.6 ; 1:bond26:1
14 19 1 0.26 23225.8 ; 1:bond27:1
19 20 1 0.232 42521 ; 1:bond28:1
20 21 1 0.344 38698 ; 1:bond29:1
21 22 1 0.48 5477.06 ; 1:bond30:1
22 23 1 0.482 4976.2 ; 1:bond31:1
23 24 1 0.482 5017.98 ; 1:bond32:1

[ angles ]
14 19 20 1 143 30 ; 1:angle1:1
19 20 21 1 109.46 53 ; 1:angle2:1
20 21 22 1 114.61 16.22 ; 1:angle3:1
21 22 23 1 143.84 19.9 ; 1:angle4:1
22 23 24 1 148.42 14.62 ; 1:angle5:1
6 15 16 1 57.87 2928 ; 1:angle6:1
5 6 15 1 50.42 3476 ; 1:angle7:1
16 17 6 1 68.77 652 ; 1:angle8:1
6 14 19 1 70 523.4 ; 1:angle9:1
3 1 2 1 84.24 500 ; 1:angle10:1
2 1 5 1 73.35 1071.6 ; 1:angle11:1
5 1 4 1 77.65 972 ; 1:angle12:1
4 1 3 1 85 400 ; 1:angle13:1
1 3 10 1 98.85 1206.22 ; 1:angle14:1
1 4 10 1 97.42 1155.18 ; 1:angle15:1
1 2 6 1 103.15 1208.72 ; 1:angle16:1
1 5 6 1 100.86 1037.2 ; 1:angle17:1
1 3 13 1 103.72 1108.8 ; 1:angle18:1
1 4 7 1 105.62 909 ; 1:angle19:1
7 5 1 1 106.01 773.28 ; 1:angle20:1
1 2 13 1 100.29 1069.6 ; 1:angle21:1
13 14 19 1 137 100 ; 1:angle22:1 mod8
6 17 18 1 143 20 ; 1:angle23:1 add mod8
18 17 16 1 114 36 ; 1:angle24:1 add mod8
4 7 8 1 47 3400 ; 1:angle25:1 add mod8
7 8 10 1 90.54 1900 ; 1:angle26:1 add mod8
6 14 13 1 91 1900 ; 1:angle27:1 add mod8
7 8 9 1 57.3 226 ; 1:angle28:1 add mod8
10 11 13 1 91 2000 ; 1:angle29:1 add mod8
10 11 12 1 65 398 ; 1:angle30:1 add mod8
13 11 12 1 148 152 ; 1:angle31:1 add mod8

[ dihedrals ]
17 16 15 7 2 51 3 ; 1:dihedral1:1
16 15 7 8 2 -57 1.34 ; 1:dihedral2:1
1 2 5 4 2 28.65 100 ; 1:dihedral3:1
1 5 4 3 2 45.78 80 ; 1:dihedral4:1
1 4 3 2 2 33.81 145 ; 1:dihedral5:1
1 3 2 5 2 43.27 150 ; 1:dihedral6:1
2 6 1 13 2 8 150 ; 1:dihedral7:1
4 10 1 7 2 12 90 ; 1:dihedral8:1
5 7 1 6 2 -5.73 50 ; 1:dihedral9:1
3 13 1 10 2 2.86 300 ; 1:dihedral10:1
10 3 1 2 2 171.35 267.52 ; 1:dihedral11:1
3 1 2 6 2 -146.13 288.06 ; 1:dihedral12:1
10 4 1 5 2 162.75 126.46 ; 1:dihedral13:1
4 1 5 6 2 171.35 167.52 ; 1:dihedral14:1
13 2 1 5 2 123.21 223.36 ; 1:dihedral15:1
2 1 5 7 2 -147.85 190.66 ; 1:dihedral16:1
13 3 1 4 2 -143.27 161.36 ; 1:dihedral17:1
3 1 4 7 2 110 80 ; 1:dihedral18:1
10 3 13 2 2 -165.61 85.58 ; 1:dihedral19:1
13 2 6 5 2 -97.42 153.18 ; 1:dihedral20:1
6 5 7 4 2 -128.94 75.08 ; 1:dihedral21:1
7 4 10 3 2 -114.6 90.64 ; 1:dihedral22:1
1 4 7 8 2 164.47 327.8 ; 1:dihedral23:1
1 5 6 15 2 163.9 266.4 ; 1:dihedral24:1
1 2 13 14 2 169.05 208.82 ; 1:dihedral25:1
1 3 10 11 2 180 157.1 ; 1:dihedral26:1
6 17 15 5 2 -28.65 24 ; 1:dihedral27:1
5 7 4 8 2 175 160 ; 1:dihedral28:1
5 6 15 16 2 172 188 ; 1:dihedral29:1
4 7 8 9 2 167.9 60 ; 1:dihedral30:1



after addition of CLA i got error. I did not found CLA CG gro file at martini website. i asked this CLA gro and itp files from my friend.

Please suggest me something. i am still stuck at this error.

Please Log in or Create an account to join the conversation.

More
9 months 2 weeks ago #8851 by vainikka
I'm not surprised, the topology looks weird to say the least.

You can find a Martini 2 compatible topology of CLA and other PSII cofactors here: cgmartini.nl/images/parameters/ITP/topologies.zip

Please Log in or Create an account to join the conversation.

More
9 months 3 days ago #8863 by saini
Thank you so much.
I solved this error.

Please Log in or Create an account to join the conversation.

More
9 months 3 days ago #8864 by saini
Dear Vainikka
I am trying to insert protein-complex in lipid bilayer by uding insane script as you suggested earlier.
used command is
python insane.py -f lhcii_em.gro -o system.gro -p system.top -x 8.5 -y 8.5 -z 8 -pbc cubic -l JPPG:5 -u JPPG:5 -center -sol W

i did energy minimzation succefully then i am tring to NVT equilibraition but after 20 ps i got the following error

1751 1755 90.0 1.1055 2.7104 0.8000
2087 2091 90.0 0.8000 1.2391 0.8000
Wrote pdb files with previous and current coordinates

Step 29315, time 58.63 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 5419197.500000, max 87309664.000000 (between atoms 170 and 171)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
159 160 55.1 0.2476 0.2457 0.2400
159 161 90.5 0.9260 0.8441 0.2400
160 161 90.4 0.9654 0.9391 0.2400
184 185 89.5 0.2397 0.4647 0.2400
184 186 89.8 1.2760 2.6281 0.2400
185 186 90.0 1.3866 2.8073 0.2400
195 196 87.3 5.6650 3.4856 0.2500
199 200 90.3 1.8829 1.6093 0.2400
199 201 91.8 1.8599 0.2164 0.2400

i am using the martini_nvt.mdp is

; VARIOUS PREPROCESSING OPTIONS =
title = Martini
;define = -DPOSRES
; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit = 0.0
dt = 0.002
nsteps = 25000 ;= 50ps
; number of steps for center of mass motion removal =
nstcomm = 100
comm-grps = Protein_lipids sol

; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 0
nstvout = 0
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 5000
nstenergy = 5000
; Output frequency and precision for xtc file =
nstxtcout = 5000
xtc_precision = 5000
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
xtc-grps = Protein_lipids sol
; Selection of energy groups =
energygrps = Protein_lipids sol

; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist = 10
; ns algorithm (simple or grid) =
ns_type = grid
; Periodic boundary conditions: xyz or none =
pbc = xyz
; nblist cut-off =
rlist = 1.2

; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = Cut-off
rcoulomb_switch = 0.0
rcoulomb = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r = 15
; Method for doing Van der Waals =
vdw_type = Cut-off
cutoff-scheme = verlet
; cut-off lengths =
rvdw_switch = 0.9
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = v-rescale ;Berendsen
; Groups to couple separately =
tc-grps = Protein_lipids sol
; Time constant (ps) and reference temperature (K) =
tau_t = 2.0 2.0
ref_t = 293 293

; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel = yes
gen_temp = 293
gen_seed = 548628

; OPTIONS FOR BONDS =
constraints = none
constraint_algorithm = Lincs
; Highest order in the expansion of the constraint coupling matrix =
lincs_order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs_warnangle = 30


please suggest me something. is their any other method available for protein insertion in bilayer?

Please Log in or Create an account to join the conversation.

More
9 months 2 days ago #8865 by vainikka
A few things you can try:

1. During initial steps of equilibration, replace constraints by bonds. Depending on the .itp file you have, you might be able to do this by defining something like -DFLEXIBLE in the .mdp file.

2. Run insane again in hopes you get a more relaxed starting configuration.

3. What is the minimum value of force you get from EM? If it is high, run it again.

4. Place your protein in an empty box and run energy minimization on it before giving it to insane.

5. Lower the timestep, and possibly lower tau-t values.

Please Log in or Create an account to join the conversation.

Time to create page: 0.148 seconds