CMM 10/25/2025
VERY ROUGH how-to steps -- This is just to familiarize yourself with gromacs/martini2 and
you will HAVE TO PUT THOUGHT INTO ACTUAL PRODUCTION SIMULATIONS!
Roughly follows Lipid Bilayer tutorial at
https://cgmartini.nl/docs/tutorials/Martini3/LipidsI/
Tarball for this how-to
You must have two key files in order to start this process. Assuming the molecule is DTAB: topology (martini_v2.0_DTAB_01.top):
#include "martini_v2.2.itp"
#include "martini_v2.0_DTAB_01.itp"
[ system ]
DTAB BILAYER SELF-ASSEMBLY
[ molecules ]
DTAB 128
..and an ITP file:
;
; Molecular topology and mapping of indices:
; +NC3-C1A-C2A-C3A-C4A
;
; 1 - 2 - 3 - 4 - 5
;
; Name
; n-dodecyltrimethylammonium bromide (DTAB)
;
; Alias
; DTAB
;
; Force field
; martini2
;
; Version
; 1.0
;
; Category
; Lipids
;
[moleculetype]
; molname nrexcl
DTAB 1
[atoms]
; id type resnr residu atom cgnr charge (mass)
1 Q0 1 DTAB NC3 1 1.0
2 SC1 1 DTAB C1A 2 0
3 SC1 1 DTAB C2A 3 0
4 C1 1 DTAB C3A 4 0
5 C1 1 DTAB C4A 5 0
[bonds]
; i j funct length force c.
; taken from Badhe https://doi.org/10.1007/s00894-024-05952-x
;
1 2 1 0.470 1250
2 3 1 0.470 1250
3 4 1 0.470 1250
4 5 1 0.470 1250
[ angles ]
; i j k funct angle force.c.
3 4 5 2 180.0 25.0
1 3 4 2 180.0 25.0
2 3 4 2 180.0 25.0
The values/parameters from the .itp file may be found in the literature or other files — you may use an existing topology from the Martini Lipodome as a guide. The bond/force constant/angle values for this particular .itp file were found/inferred from Badhe, et al https://doi.org/10.1007/s00894-024-05952-x.
You must decide how to translate the molecule into martini2 beads; In this case it is easy as a Q0 bead is charged, and the NC3 name is commonly used for trimethylammonium ion. C1 beads represent 4 carbons, and SC1 represent 3, so C14 = 2xSC1 + 2xC1, Even with this mapping, different amphiphiles may not turn out to be unique; see the discussion in the tutorial, and realize that Martini3 makes this much less of a problem.
Using this particular set of files as a guide, it should be relatively simple to construct the corresponding gemini surfactant.
you must have the following files to make this work:
DTAB-em.gro -- the coordinates in .gro format of your surfactant
martini_v2.0_DTAB_01.itp -- the itp file for that surfactant
martini_v2.2.itp -- the main martini2 itp file
martini_v2.0_ions.itp
martini_v2.0_Br.itp -- the itp for bromide (or you can put it into the above ions file)
martini_v2.0_solvents.itp
waterbox.gro -- equilibrated martini water box
various .mdp files as supplied in the tutorial
DTAB.top -- the topology file (easy to make)
1) assemble the DTAB. adjust box size and number depending on what system (standard disclaimer)
gmx insert-molecules -ci DTAB-em.gro -box 7.5 7.5 7.5 -nmol 128 -radius 0.21 -try 500 -o 128_noW.gro
NOTE: If you are adding polyatomic anions such as benzoate or sulfate, you must use insert-molecules successively (instead of genion, below):
gmx insert-molecules -ci DTAB-em.gro -box 10 10 10 -nmol 128 -radius 0.21 -try 500 -o 128_DTAB.gro
gmx insert-molecules -ci BzO-0.gro -f 128_DTAB.gro -box 10 10 10 -nmol 128 -radius 0.21 -try 500 -o 128_DTAB_BzO.gro
2) Add martini water using solvate:
gmx solvate -cp 128_noW.gro -cs water.gro -o waterbox.gro -maxsol 768 -radius 0.21
(because of the initial patchy association during minimization, probably the 'maxsol' should be increased to some ridiculously large number? Try this!)
50 DTAB, asked for maxsol 10000, got 7695.
50 DTAB in 100 L box; asked for 8,000,000 got 8,000,000 (!)
10 DTAB in 58.5 L box; asked for 1,600,000 got
50 DTAB in 46.4 box is 1.5 mM is 800,000 W got 795395
You could minimize at this time -- see the tutorial page. However you may (should?) wait until the bromides are added. Anytime you do actual grompp and/or mdrun, you DTAB.top will be modified to add any new species.
So, at this point, you need to add a line into DTAB.top to account for the newly-added waters ("W line only"): You need to determine what this number is by reading the previous outputs.
[ molecules ]
DTAB 128
W 640
3) add the bromide ions using genion. You could add positive or negative ions, but in this case we want to neutralize with bromides only. A good example of genion use may be found at http://www.mdtutorials.com/gmx/lysozyme/04_ions.html
You must have an "ions.mdp" despite the fact that you aren't going to actually run anything at this time.
'genion' will ask which group to affect; you should pick the W group (which is probably group 3).
gmx grompp -f ions.mdp -c waterbox.gro -p DTAB.top -o ions.tpr
gmx genion -s ions.tpr -o waterbox_ions.gro -p DTAB.top -nname BR- -neutral
If you examine your .top file, you should now see the BR- ions added to the end.
4) now minimize this DTAB + water + bromide box as specified in the Bilayers tutorial. You may want to make the output file (minimized.gro) something more descriptive. If you are running on Apple Silicon, gromacs will try and grab all the processors by default. This actually requires some thought; for now, make sure and throw the '-ntmpi 1' switch to use only 1 processor.
4A) At this point, you will probably have to (re-)create an 'index.ndx' file if you get warnings about BR- in the next step.(???) ACTUALLY this is because BR-'s group is "ION"; the .mdp file should refer to "ION" not "BR-".
This will obviate the "index.ndx" usage.
DON'T THINK THIS IS TRUE BUT SEE LATER STEP w/ combining groups (martiniglass)
You may have to manually edit the .gro file to change any mention of "BR" to "BR-". I am not sure why this happens, but it may be due to the "-" character. You may be tempted to change the '-maxwarn' switch to something greater than 1 as the processor suggests, but the solution is to make the replacements.
gmx grompp -f minimization.mdp -c waterbox_ions.gro -p DTAB.top -o DTAB-min-solvent.tpr
gmx mdrun -ntmpi 1 -s DTAB-min-solvent.tpr -v -c minimized.gro
gmx make_ndx -f minimized.gro -o DTAB_50-md.ndx
5) Now you should be able to run dynamics! (Please note you probably need to utilize the index.ndx file. NOT TRUE if system was set up properly/consistiently)
gmx grompp -f martini_md.mdp -c minimized.gro -n index.ndx -p DTAB.top -o DTAB-md.tpr
gmx mdrun -ntmpi 1 -s DTAB-md.tpr -v -x DTAB-md.xtc -c DTAB-md.gro -n index.ndx
w/ 50 DTAB, 7695 - 50 W, 50 BR- "Step 25 Warning: pressure scaling more than 1%..."
Seems OK -- micelles form after 100000 steps 3 ns, so probably (well?) above CMC
You should read the rest of the tutorial -- visualization is a whole 'nother can of worms. Using MartiniGlass can work well. Look for a how-to on this topic.
Visualization with martiniglass
Make sure it is installed via pip.
https://martiniglass.readthedocs.io/en/latest/tutorials/non_protein.html
https://martiniglass.readthedocs.io/en/latest/advanced_options.html#waterless-trajectories
should remove water
need to use make_ndx to gather the non-water residues into one group:
gmx make_ndx -f DTAB_50-md.gro -o DTAB_50-md.ndx
Inside the "editor" choose
> n1 | n2
(this will join groups n1 and n2 into a new group -- choose the DTAB and ION groups, for example)
Now run martiniglass. We will use "vis.xxx" as the filenames, but of course you may choose anything:
martiniglass -p DTAB_mic_50.top -f DTAB_50-md.gro -vf
gmx trjconv -f DTAB_50-md.gro -s DTAB_50-md.gro -n DTAB_50-md.ndx -o vis.gro
You will be asked which groups; you should choose the combined group (called DTAB_ION)
At this time you may also want to create a trajectory file with the same residues as these files:
gmx trjconv -f DTAB_50-md.xtc -s DTAB_50-md.gro -n DTAB_50-md.ndx -o vis.xtc
Now, assuming you have vmd aliased in your shell, you run it:
vmd vis.gro -e vis.vmd
Martiniglass by default creates many representations; most are uneeded. Find the representation with the lipids (POPC, POPE, etc) and change all those names to "DTAB". You should see the polar tetraammonium head as purple, and the HC tail(s) as grey caterpillars.
There is a problem with visualizing trajectories. Sometimes atoms (beads) cross the periodic boundaries of the box and the molecule appears to stretch across the whole box, giving the appearance of box-length rds. This can be fixed by a form of periodic recentering (the proper -pbc setting was found by trial-and-error):
gmx trjconv -f DTAB_50-md.xtc -s DTAB_50-md.tpr -n DTAB_50-md.ndx -pbc mol -ur compact -o vis-compact.xtc
(If you have problems of this sort, it is a good idea to test different combinations of the PBC/centering commands (of course .trr could be .xtc or .gro)
Commands are as follows:
gmx trjconv -s md.tpr -f md.trr -o md1.trr -center -pbc whole -n index.ndx
gmx trjconv -s md.tpr -f md1.trr -o md2.trr -pbc nojump -n index.ndx
gmx trjconv -s md.tpr -f md2.trr -o md3.trr -pbc mol -ur compact -n index.ndx
This trajectory file can be read into the previous (static) vmd molecule (which utilized the single-frame .gro file) in order to produce a trajectory and/or movies.
VERY ROUGH how-to steps -- This is just to familiarize yourself with gromacs/martini2 and
you will HAVE TO PUT THOUGHT INTO ACTUAL PRODUCTION SIMULATIONS!
Roughly follows Lipid Bilayer tutorial at
https://cgmartini.nl/docs/tutorials/Martini3/LipidsI/
Tarball for this how-to
PRELUDE: Creating the primary files
You must have two key files in order to start this process. Assuming the molecule is DTAB: topology (martini_v2.0_DTAB_01.top):
#include "martini_v2.2.itp"
#include "martini_v2.0_DTAB_01.itp"
[ system ]
DTAB BILAYER SELF-ASSEMBLY
[ molecules ]
DTAB 128
..and an ITP file:
;
; Molecular topology and mapping of indices:
; +NC3-C1A-C2A-C3A-C4A
;
; 1 - 2 - 3 - 4 - 5
;
; Name
; n-dodecyltrimethylammonium bromide (DTAB)
;
; Alias
; DTAB
;
; Force field
; martini2
;
; Version
; 1.0
;
; Category
; Lipids
;
[moleculetype]
; molname nrexcl
DTAB 1
[atoms]
; id type resnr residu atom cgnr charge (mass)
1 Q0 1 DTAB NC3 1 1.0
2 SC1 1 DTAB C1A 2 0
3 SC1 1 DTAB C2A 3 0
4 C1 1 DTAB C3A 4 0
5 C1 1 DTAB C4A 5 0
[bonds]
; i j funct length force c.
; taken from Badhe https://doi.org/10.1007/s00894-024-05952-x
;
1 2 1 0.470 1250
2 3 1 0.470 1250
3 4 1 0.470 1250
4 5 1 0.470 1250
[ angles ]
; i j k funct angle force.c.
3 4 5 2 180.0 25.0
1 3 4 2 180.0 25.0
2 3 4 2 180.0 25.0
The values/parameters from the .itp file may be found in the literature or other files — you may use an existing topology from the Martini Lipodome as a guide. The bond/force constant/angle values for this particular .itp file were found/inferred from Badhe, et al https://doi.org/10.1007/s00894-024-05952-x.
You must decide how to translate the molecule into martini2 beads; In this case it is easy as a Q0 bead is charged, and the NC3 name is commonly used for trimethylammonium ion. C1 beads represent 4 carbons, and SC1 represent 3, so C14 = 2xSC1 + 2xC1, Even with this mapping, different amphiphiles may not turn out to be unique; see the discussion in the tutorial, and realize that Martini3 makes this much less of a problem.
Using this particular set of files as a guide, it should be relatively simple to construct the corresponding gemini surfactant.
Building the simulation box
you must have the following files to make this work:
DTAB-em.gro -- the coordinates in .gro format of your surfactant
martini_v2.0_DTAB_01.itp -- the itp file for that surfactant
martini_v2.2.itp -- the main martini2 itp file
martini_v2.0_ions.itp
martini_v2.0_Br.itp -- the itp for bromide (or you can put it into the above ions file)
martini_v2.0_solvents.itp
waterbox.gro -- equilibrated martini water box
various .mdp files as supplied in the tutorial
DTAB.top -- the topology file (easy to make)
1) assemble the DTAB. adjust box size and number depending on what system (standard disclaimer)
gmx insert-molecules -ci DTAB-em.gro -box 7.5 7.5 7.5 -nmol 128 -radius 0.21 -try 500 -o 128_noW.gro
NOTE: If you are adding polyatomic anions such as benzoate or sulfate, you must use insert-molecules successively (instead of genion, below):
gmx insert-molecules -ci DTAB-em.gro -box 10 10 10 -nmol 128 -radius 0.21 -try 500 -o 128_DTAB.gro
gmx insert-molecules -ci BzO-0.gro -f 128_DTAB.gro -box 10 10 10 -nmol 128 -radius 0.21 -try 500 -o 128_DTAB_BzO.gro
Solvating with water beads
2) Add martini water using solvate:
gmx solvate -cp 128_noW.gro -cs water.gro -o waterbox.gro -maxsol 768 -radius 0.21
(because of the initial patchy association during minimization, probably the 'maxsol' should be increased to some ridiculously large number? Try this!)
50 DTAB, asked for maxsol 10000, got 7695.
50 DTAB in 100 L box; asked for 8,000,000 got 8,000,000 (!)
10 DTAB in 58.5 L box; asked for 1,600,000 got
50 DTAB in 46.4 box is 1.5 mM is 800,000 W got 795395
You could minimize at this time -- see the tutorial page. However you may (should?) wait until the bromides are added. Anytime you do actual grompp and/or mdrun, you DTAB.top will be modified to add any new species.
So, at this point, you need to add a line into DTAB.top to account for the newly-added waters ("W line only"): You need to determine what this number is by reading the previous outputs.
[ molecules ]
DTAB 128
W 640
Adding bromide counter-ions
3) add the bromide ions using genion. You could add positive or negative ions, but in this case we want to neutralize with bromides only. A good example of genion use may be found at http://www.mdtutorials.com/gmx/lysozyme/04_ions.html
You must have an "ions.mdp" despite the fact that you aren't going to actually run anything at this time.
'genion' will ask which group to affect; you should pick the W group (which is probably group 3).
gmx grompp -f ions.mdp -c waterbox.gro -p DTAB.top -o ions.tpr
gmx genion -s ions.tpr -o waterbox_ions.gro -p DTAB.top -nname BR- -neutral
If you examine your .top file, you should now see the BR- ions added to the end.
Minimization
4) now minimize this DTAB + water + bromide box as specified in the Bilayers tutorial. You may want to make the output file (minimized.gro) something more descriptive. If you are running on Apple Silicon, gromacs will try and grab all the processors by default. This actually requires some thought; for now, make sure and throw the '-ntmpi 1' switch to use only 1 processor.
4A) At this point, you will probably have to (re-)create an 'index.ndx' file if you get warnings about BR- in the next step.(???) ACTUALLY this is because BR-'s group is "ION"; the .mdp file should refer to "ION" not "BR-".
This will obviate the "index.ndx" usage.
DON'T THINK THIS IS TRUE BUT SEE LATER STEP w/ combining groups (martiniglass)
You may have to manually edit the .gro file to change any mention of "BR" to "BR-". I am not sure why this happens, but it may be due to the "-" character. You may be tempted to change the '-maxwarn' switch to something greater than 1 as the processor suggests, but the solution is to make the replacements.
gmx grompp -f minimization.mdp -c waterbox_ions.gro -p DTAB.top -o DTAB-min-solvent.tpr
gmx mdrun -ntmpi 1 -s DTAB-min-solvent.tpr -v -c minimized.gro
gmx make_ndx -f minimized.gro -o DTAB_50-md.ndx
Running dynamics
5) Now you should be able to run dynamics! (Please note you probably need to utilize the index.ndx file. NOT TRUE if system was set up properly/consistiently)
gmx grompp -f martini_md.mdp -c minimized.gro -n index.ndx -p DTAB.top -o DTAB-md.tpr
gmx mdrun -ntmpi 1 -s DTAB-md.tpr -v -x DTAB-md.xtc -c DTAB-md.gro -n index.ndx
w/ 50 DTAB, 7695 - 50 W, 50 BR- "Step 25 Warning: pressure scaling more than 1%..."
Seems OK -- micelles form after 100000 steps 3 ns, so probably (well?) above CMC
You should read the rest of the tutorial -- visualization is a whole 'nother can of worms. Using MartiniGlass can work well. Look for a how-to on this topic.
VISUALIZATION — Martiniglass and VMD
Visualization with martiniglass
Make sure it is installed via pip.
https://martiniglass.readthedocs.io/en/latest/tutorials/non_protein.html
https://martiniglass.readthedocs.io/en/latest/advanced_options.html#waterless-trajectories
should remove water
need to use make_ndx to gather the non-water residues into one group:
gmx make_ndx -f DTAB_50-md.gro -o DTAB_50-md.ndx
Inside the "editor" choose
> n1 | n2
(this will join groups n1 and n2 into a new group -- choose the DTAB and ION groups, for example)
Now run martiniglass. We will use "vis.xxx" as the filenames, but of course you may choose anything:
martiniglass -p DTAB_mic_50.top -f DTAB_50-md.gro -vf
gmx trjconv -f DTAB_50-md.gro -s DTAB_50-md.gro -n DTAB_50-md.ndx -o vis.gro
You will be asked which groups; you should choose the combined group (called DTAB_ION)
At this time you may also want to create a trajectory file with the same residues as these files:
gmx trjconv -f DTAB_50-md.xtc -s DTAB_50-md.gro -n DTAB_50-md.ndx -o vis.xtc
Now, assuming you have vmd aliased in your shell, you run it:
vmd vis.gro -e vis.vmd
Martiniglass by default creates many representations; most are uneeded. Find the representation with the lipids (POPC, POPE, etc) and change all those names to "DTAB". You should see the polar tetraammonium head as purple, and the HC tail(s) as grey caterpillars.
There is a problem with visualizing trajectories. Sometimes atoms (beads) cross the periodic boundaries of the box and the molecule appears to stretch across the whole box, giving the appearance of box-length rds. This can be fixed by a form of periodic recentering (the proper -pbc setting was found by trial-and-error):
gmx trjconv -f DTAB_50-md.xtc -s DTAB_50-md.tpr -n DTAB_50-md.ndx -pbc mol -ur compact -o vis-compact.xtc
(If you have problems of this sort, it is a good idea to test different combinations of the PBC/centering commands (of course .trr could be .xtc or .gro)
Commands are as follows:
gmx trjconv -s md.tpr -f md.trr -o md1.trr -center -pbc whole -n index.ndx
gmx trjconv -s md.tpr -f md1.trr -o md2.trr -pbc nojump -n index.ndx
gmx trjconv -s md.tpr -f md2.trr -o md3.trr -pbc mol -ur compact -n index.ndx
This trajectory file can be read into the previous (static) vmd molecule (which utilized the single-frame .gro file) in order to produce a trajectory and/or movies.