Solvate an Initial Configuration

Classes

class solvation_equilibration.System(build_monomer, weight_percent, ratio, tolerance=1, solute='HOH', nopores=4, ncolumns=5, monomers_per_column=20, p2p=4.5, parallel_displaced=0, dbwl=0.37, random_seed=0, mpi=False, nproc=4)

__init__(build_monomer, weight_percent, ratio, tolerance=1, solute='HOH', nopores=4, ncolumns=5, monomers_per_column=20, p2p=4.5, parallel_displaced=0, dbwl=0.37, random_seed=0, mpi=False, nproc=4)

Get unit cell build parameters and determine the number of water molecules required for each region in order to have the correct final composition.

Parameters:

• build_monomer (str) – Name of liquid crystal monomer to use to build unit cell (no file extension)
• weight_percent (float) – Weight percent of water in the whole system
• ratio (float) – Ratio of water in the pores to water in the tails
• tolerance (int) – Acceptable error in the total number of water molecules placed in the pore region
• solute (str) – name of solute used to solvate system
• nopores (int) – Number of pores in unit cell
• ncolumns (int) – Number of stacked monomer columns per pore
• monomers_per_column (int) – Number of monomers stacked into a column
• p2p (float) – Distance between pores (nm)
• parallel_displaced (float) – Angle of wedge formed between line extending from pore center to monomer and line from pore center to vertically adjacent monomer head group
• dbwl (float) – Distance between stacked monomers (nm)
• random_seed (int) – Monomer columns are randomly displaced in the z-direction when the initial configuration is built. Specify a random seed so that the same displacement is achieved each time the radius is changed. I think this helps with convergence, but it hasn’t been extensively tested.
• mpi (bool) – Run the MPI version of GROMACS
• np (int) – number of MPI process if mpi = True

calculate_pore_water(config)

Determine the total number of water molecules within the pore radius. Update database with system configuration parameters.

Parameters:config (str) – Name of .gro configuration of which to calculate pore water content

equilibrate(input_files=True, length=50, force=1000000, restraint_residue='HII', restraint_axis='xyz', ring_restraints=('C', 'C1', 'C2', 'C3', 'C4', 'C5'))

Simulate the unit cell with restraints placed on the head group

Parameters:

• input_files (bool) – Generate GROMACS .mdp and topology files
• length (int) – Simulation length (ps)
• force (float or int) – Force with which to restrain head groups kJ mol^-1 nm^-2
• restraint_residue (str) – Name of residue to which position restraints will be applied
• restraint_axis (str) – Axis/axes along which head groups should be restrained
• ring_restraints (tuple) – Names of head group atoms

full_equilibration(forces, fully_solvated='solvated_final.gro', l_nvt=50, l_berendsen=5000, l_pr=400000, restraint_residue='HII', restraint_axis='xyz', ring_restraints=('C', 'C1', 'C2', 'C3', 'C4', 'C5'))

Simulate the unit cell with a sequence of decreasing restraints placed on the head group

Parameters:

• forces (tuple or list) – sequence of forces to apply to ring_restraints (\(\dfrac{kJ}{mol~nm^2}\))
• fully_solvated (str) – name of fully solvated coordinate file
• l_nvt (int) – Length of short restrained simulations (ps)
• l_berendsen (int) – Length of equilibration simulation run with berendsen pressure control (ps)
• l_pr (int) – Length of long equilibration simulation run with Parrinello-Rahman pressure control (ps)
• restraint_residue (str) – Name of residue to which position restraints will be applied
• restraint_axis (str) – Axis/axes along which head groups should be restrained
• ring_restraints (tuple) – Names of head group atoms

place_water_tails(output)

Place the desired number of water molecules in the tails. This is done by randomly inserting water molecules in close proximity to the tails, one-by-one. A short energy minimzation is performed between each insertion.

Parameters:output (str) – Name of final configuration

query_database(database='water_content.db', guess_range=(0.4, 1), guess_stride=0.2)

Read an SQL database, of number of water molecules versus pore radius, to inform the next choice of pore radius. The decision process bisects the two radii with water contents closest to the desired concentration. If there is a single data point, an educated guess is made. If there are no data points, then make a random guess.

Parameters:

• database (str) – Name of database
• guess_range (tuple) – If water_content.db has no entries for the build monomer, an initial radius will be randomly selected from this range
• guess_stride (float) – How far above/below the highest/lowest value to make the next guess at pore radius if you need more/less water than the bounds of the water

write_final_pore_configuration()

Write the configuration with the correct number of water molecules placed in the pore region, then rwrite topology files.

Examples

# This pseudocode illustrates the commands that will build and solvate a system from scratch
sys = System()

while not sys.converged:
    sys.query_database()  # make a guess at the pore radius
    sys.equilibrate()  # Build system with guess radius, add water and make sure it is stable
    sys.calculate_pore_water() # Determine the total water in the pore. Stop loop if within tolerance

sys.write_final_pore_configuration()  # write out coordinate file with correct no. of pore water
sys.place_water_tails()  # put the appropriate amount of water in the tails
sys.full_equilibration()  # Fully equilibrate the fully solvated system

Command Line Functionality

For convenience, a full equilibration can be achieved using a single python script which contains the above class and runs the sequence of commands in the example code block above.

Solvate system and adjust so there is a certain wt % of water

usage: solvation_equilibration.py [-h] [-ratio RATIO] [-wt WEIGHT_PERCENT]
                                  [-tol TOLERANCE]
                                  [-guess_range GUESS_RANGE [GUESS_RANGE ...]]
                                  [-guess_stride GUESS_STRIDE] [-o OUTPUT]
                                  [-seed RANDOM_SEED] [-mpi] [-np NPROC]
                                  [-b BUILD_MONOMER] [-m MONOMERS_PER_COLUMN]
                                  [-c NCOLUMNS] [-r PORE_RADIUS] [-p P2P]
                                  [-n NOPORES] [-d DBWL]
                                  [-pd PARALLEL_DISPLACED]
                                  [-angles ANGLES [ANGLES ...]]
                                  [-ring_restraints RING_RESTRAINTS [RING_RESTRAINTS ...]]
                                  [-forces FORCES]
                                  [--restraint_residue RESTRAINT_RESIDUE]
                                  [--restraint_axis RESTRAINT_AXIS]
                                  [-l_nvt LENGTH_NVT] [-lb LENGTH_BERENDSEN]
                                  [-lpr LENGTH_PARRINELLO_RAHMAN]

Named Arguments

-ratio, --ratio
	Ratio of water in pores to water in tails Default: 2
-wt, --weight_percent
	Total weight percent of water Default: 10
-tol, --tolerance
	Number of water molecules Default: 1
-guess_range	If water_content.db has no entries for the build monomer, an initial radius will be randomly selected from this range Default: [0.4, 1]
-guess_stride	How far above/below the highest/lowest value tomake the next guess at pore radius if you need more/less water than the bounds of the water content database (nm) Default: 0.2
-o, --output	Name of fully solvated output file Default: “solvated_final.gro”
-seed, --random_seed
	Numpy random seed Default: 0
-mpi, --mpi	Run MD simulations in parallel Default: False
-np, --nproc	Number of MPI processes Default: 4
-b, --build_monomer
	Name of single monomer used to build full system Default: “NAcarb11V”
-m, --monomers_per_column
	Number of monomers to stack in eachcolumn Default: 20
-c, --ncolumns	Number of columns used to build each pore Default: 5
-r, --pore_radius
	Initial guess at pore radius (nm) Default: 0.6
-p, --p2p	Initial pore-to-pore distance (nm) Default: 4.5
-n, --nopores	Number of pores (only works with 4 currently) Default: 4
-d, --dbwl	Distance between vertically stacked monomers(nm) Default: 0.37
-pd, --parallel_displaced
	Angle of wedge formed between lineextending from pore center to monomer and line from pore center to vertically adjacent monomerhead group. Default: 0
-angles, --angles
	Angles betweenbox vectors Default: [90, 90, 60]
-ring_restraints, --ring_restraints
	Name of atoms used to restrain head groups during initial equilibration. Default: [‘C’, ‘C1’, ‘C2’, ‘C3’, ‘C4’, ‘C5’]
-forces	Sequence of force constants toapply to ring restraints Default: [1000000, 3162, 56, 8, 3, 2, 1, 0]
--restraint_residue
	Name of residue to which ring_restraintatoms belong Default: “HII”
--restraint_axis
	Axes along which to apply position restraints Default: “xyz”
-l_nvt, --length_nvt
	Length of restrained NVT simulations (ps) Default: 50
-lb, --length_berendsen
	Length of simulation using berendsenpressure control Default: 5000
-lpr, --length_Parrinello_Rahman
	Length of simulation using Parrinello-Rahman pressure control Default: 400000