Relaxation to Build Stable Membrane Systems and Running MD

by Stephan Schott-Verdugo (s.schott-verdugo@fz-juelich.de)

Minimizing and Equilibrating a packed membrane system

Once you have built your membrane systems, you will want to create a stable system to simulate. Here you can find a list of input include for minimizing and equilibrating the structures that will be generated, irrespective of the membrane composition or if a protein is included. You can also use these scripts with systems packed with CHARMM-GUI.

The general minimization workflow used here consists of 5 steps (all of which I personally run in CPUs, hopefully using the MPI implementation, particularly if the packing is complex; chances of things failing in the minimization step are, from higher to lower likelihood, and from faster to slower: pmemd.cuda>pmemd.MPI>pmemd>sander). In addition, each step uses both, the steepest descent and the conjugated gradient algorithm; we have seen that you need to use both to remove clashes in most cases!:

01_min_water_primary.in       Initial water minimization, restraining the membrane residues (and protein residues if present) with 25 kcal mol^-1Å^-2
02_min_water_secondary.in   Water minimization, restraining the membrane residues (and protein residues if present) with 5 kcal mol^-1Å^-2
03_min_primary.in                  Membrane minimization, (restraining the protein and ligand residues, if present, with 5 kcal mol^-1Å^-2)
04_min_secondary.in              Membrane minimization, (restraining the protein residues, if present, with 5 kcal mol^-1Å^-2)
05_min_all.in                          Minimization without restraints

In all these include, you have to make sure to replace the restraining mask for the proper residue span in the system you generated. For example, if you have a system with a protein of 100 residues, and you have 70 lipids per leaflet, you will have to replace :PROTEIN+LIPID for the corresponding mask that selects all of them (considering that each lipid has 3 "residues", and we have 2 leaflets, the selection will probably be :1-520). The previous minimization might seem long, but in practice, it ensures to remove as many clashes as possible, and if the system reaches a minimum, it will exit the minimization procedure with a LINMIN failure (completely normal!). In case of failure (positive energies, overlapping hydrogens, etc), you can as well check the dx0 flag.

You can use the previous scripts with something like the following:

DO_PARALLEL=YOUR_PARALLEL_CMD_HERE EXE=pmemd.MPI PRMTOP=YOUR_TOP_HERE INPCRD=YOUR_CRD_HERE OLD=$INPCRD NEW=01_min_water_primary $DO_PARALLEL $EXE -O -i $NEW.in -o $NEW.out -p $PRMTOP -c $OLD -r $NEW.restrt -ref $OLD OLD=$NEW.restrt NEW=02_min_water_secondary $DO_PARALLEL $EXE -O -i $NEW.in -o $NEW.out -p $PRMTOP -c $OLD -r $NEW.restrt -ref $OLD OLD=$NEW.restrt NEW=03_min_lipids_primary $DO_PARALLEL $EXE -O -i $NEW.in -o $NEW.out -p $PRMTOP -c $OLD -r $NEW.restrt -ref $OLD OLD=$NEW.restrt NEW=04_min_lipids_secondary $DO_PARALLEL $EXE -O -i $NEW.in -o $NEW.out -p $PRMTOP -c $OLD -r $NEW.restrt -ref $OLD OLD=$NEW.restrt NEW=05_min_all $DO_PARALLEL $EXE -O -i $NEW.in -o $NEW.out -p $PRMTOP -c $OLD -r $NEW.restrt -ref $OLD

Once you have a minimized your system, we have to thermalize it, and let the density reach a desirable value. Running these steps in GPUs should cause no issue, as long as all clashes were removed in the previous steps. By default, PACKMOL-Memgen uses the LEaP setBox vdw method to set the box dimensions, which can leave a considerable space between periodic images, and hence a very low initial density (Note: you might prefer to use the --tight_box flag if you want to force box dimensions according to the estimated area per lipid, though this might cause unresolved clashes in the periodic boundary during the minimization). This is rapidly solved as long as we do not simulate too long using the NVT ensemble in the initial steps:

06_nvt_initial.in    Very short (5 ps) NVT heating from 0 to 100K (with restraints on the protein, and optionally on the membrane). This primes the system velocities.
07_npt_initial.in    Initial (115 ps) NPT heating from 100 to 300K (with restraints on the protein). This will let the box dimensions adapt to the system shape.
08_npt_pbceq.in    PBC equilibration, completing 5 ns of simulation in NPT.

As for the minimization, you can run the thermalization step using something like the following:

DO_PARALLEL=YOUR_PARALLEL_CMD_HERE EXE=pmemd.cuda PRMTOP=YOUR_TOP_HERE INPCRD=05_min_all.restrt OLD=$INPCRD NEW=06_nvt_initial $DO_PARALLEL $EXE -O -i $NEW.in -o $NEW.out -p $PRMTOP -c $OLD -r $NEW.restrt -ref $OLD -x $NEW.nc OLD=$NEW.restrt NEW=07_ntp_initial $DO_PARALLEL $EXE -O -i $NEW.in -o $NEW.out -p $PRMTOP -c $OLD -r $NEW.restrt -ref $OLD -x $NEW.nc OLD=$NEW.restrt NEW=08_ntp_pbceq $DO_PARALLEL $EXE -O -i $NEW.in -o $NEW.out -p $PRMTOP -c $OLD -r $NEW.restrt -ref $OLD -x $NEW.nc

Please be aware that to properly equilibrate the system you might require longer than the 5 ns performed here; membranes seem to stabilize their thickness and area per lipid values at > 50 ns!