Modifying topol.top with each step of peptide insertion over membrane from CHARMM-GUI

User discussions
,
1

GROMACS version: 2025.2
GROMACS modification: No

I got membrane from charmm-gui and equilibrated it for 5 ns with 4 steps of NVT and 2 steps of NPT. No matter what I do the waters above membrane always shrink. So I decided to remove waters from the equilibrated membrane, added some peptides and resolvated the system.

For peptide,

gmx pdb2gmx -f peptide.pdb -o peptide.gro -p peptide.top -i peptide.itp
Select charmm36 and TIP3P , same as the membrane system in Charmm-GUI.

Then, I ran
gmx edifconf -f membrane.gro -o file.gro -c -box x y z (last line of membrane.gro file)
gmx make_ndx -f file.gro -n index.ndx ;Selected the entire system
gmx -trjconv -f file.gro -o new.gro -pbc mol -n index.ndx ; (selected the system) To keep the system at bottom as I expand the box only above the membrane.

Then,
gmx make_ndx -f new.gro -n index.ndx ; Selected only the membrane systems

; Only keeps the membrane in the box, removes everything else)
gmx editconf -f new.gro -o only_membrane.gro -n index.ndx
I deleted the ions and water from my topol.top file

: Resolvate
gmx solvate -cp only_membrane.gro -cs spc216.gro -o solvated.gro
topol.top changed as usual.

; Adding peptides

gmx insert-molecules -f out.gro -ci peptide.gro -o pept_mem.gro -nmol 25 -replace -rot z -ip positions.dat
Only 11 molecules were inserted.

Things get complicated here. So I add the peptide numbers to topol.top like this:

;;
;; Generated by CHARMM-GUI FF-Converter
;;
;; Correspondance:
;; jul316@lehigh.edu or wonpil@lehigh.edu
;;
;; The main GROMACS topology file
;;

; Include forcefield parameters
#include “toppar2/forcefield.itp”
#include “toppar/POPG.itp”
#include “toppar/TIP3P.itp”
#include “peptide.top”

[ system ]
; Name
Title in water

**[ molecules ]
; Compound #mols
POPG 460
DCLE 2277
SOL 53900
peptide 11
**
After that I try to minimize the system with :
gmx grompp -f mdout.mdp -c pept_mem.gro -p topol.top -o minim.tpr

But I keep getting the error:

Fatal error:
Syntax error - File forcefield.itp, line 11
Last line read:
‘[ defaults ]’
Invalid order for directive defaults

image
image856×717 201 KB

2

You cannot #include a system topology. This topology will call the force field that you specify when running pdb2gmx. So you specify toppar2/forcefield.itp and then you are also calling charmm36.ff/forcefield.itp from within peptide.top. You thus have two [defaults] directives, which are the core instructions for the energy function. This is syntactically not allowed. You have to remove all references to inclusion of a force field and any system-level directives (e.g., [system] and [molecules]) from peptide.top to make a proper .itp file that can be included.

3

Thank you, I followed your instructions and it worked out. But I am having more errors in subsequent steps.
When I use this mdp to neutralize the system I am getting some notes.
minim.mdp:
; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = cutoff ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; Short-range electrostatic cut-off
rvdw = 1.2 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions

Outputs:
Command line:
gmx grompp -f minim.mdp -o minim.tpr -c solv.gro -p topol.top -n index.ndx

Ignoring obsolete mdp entry ‘ns_type’

NOTE 1 [file minim.mdp]:
With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
that with the Verlet scheme, nstlist has no effect on the accuracy of
your simulation.

Setting the LD random seed to -281059331

Generated 167799 of the 167910 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 1

Generated 117432 of the 167910 1-4 parameter combinations

Excluding 3 bonded neighbours molecule type ‘POPG’

Excluding 3 bonded neighbours molecule type ‘DCLE’

Excluding 3 bonded neighbours molecule type ‘Protein’

Excluding 2 bonded neighbours molecule type ‘SOL’

NOTE 2 [file topol.top, line 25]:
System has non-zero total charge: -388.000000
Total charge should normally be an integer. See

for discussion on how close it should be to an integer.

Number of degrees of freedom in T-Coupling group rest is 473517.00
The integrator does not provide a ensemble temperature, there is no system ensemble temperature

NOTE 3 [file minim.mdp]:
You are using a plain Coulomb cut-off, which might produce artifacts.
You might want to consider using PME electrostatics.

But it works.

Then I run energy minimization after this with this mdp:

integrator = steep
emtol = 1000.0
nsteps = 5000
nstlist = 10
cutoff-scheme = Verlet
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = PME
rcoulomb = 1.2
;
constraints = h-bonds
constraint_algorithm = LINCS

Then I get some errors while running
gmx grompp -f ions.mdp -o neutral.tpr -c neutral.gro -p topol.top -n index.ndx
gmx mdrun -v -deffnm neutral

And these were in the results, although the system converged to the desired steepest descent.
step 18: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.

Step 34, time 0.034 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000066, max 0.010883 (between atoms 7281 and 7280)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
7282 7280 35.2 0.1110 0.1108 0.1111