Note: Best viewed in the notebook viewer
Protocol for initial configuration for confinedBSA
Meetings
Example meeting
(Top)
(Top)
IPTS-18245 - proposal document
Two samples:
Susceptibility plots:
| DBSA |
HBSA |
Comparison of MSD (label says 'Elastic Intensity' but it's wrong):
Maximum MSD corresponds to ~0.4A, very minimal flucuation. Thus, there is very little deuterated water surrounding the protein in the silica-encapsulated DBSA sample.
Slides that Panchao sent me providing more info on the samples:
The last slide shows the relative weights of BSA (66054.4 amu per monomer), silica (60.08 amu per SiO2 molecule), and water (18.015 amu per molecule), corresponding to 1545 SiO2 and 622 water molecules for each BSA monomer.
(Top)
References for mean square displacements in bulk materials containing silica
Simulations of proteins interacting with silica
(Top)
. *.top: (ASCII) The topology contains all the information necessary to define the molecule within a simulation. Thi- information includes nonbonded parameters (atom types and charges) as well as bonded parameters (bonds, angles, and dihedrals). Does not contain coordinates for the atoms.
*.gro: (ASCII) GROMACS-formatted structure file that contains residue and atom names of our system with the nomenclature of the specified force field. Also coordinates.
*.itp: (ASCII) Indicates a harmonic restraint of specified force constant to be applied for certain atoms of the simulation.
*.mdp: (ASCII) Run configuration parameters.
*.tpr: (binary) Run input. Integrates information from the .top, .gro, and *.mdp files.
*.log: (ASCII) standard output of the Gromacs run.
*.edr: (binary) values for various thermodynamic properties of the system versus time.
*.trr: (binary) trajectory (coordinates, and optionally velocities) versus time, in full-precision.
*.cpt: (binary) Checkpoint containing all the state variables, including coordinates and velocities, that allow one to continue a simulation.
*.xvg: (ASCII) Output of gmx energy, a file to be read by the plotting tool Xmgrace.
*.ndx: (ASCII) index containing user-defined sets of atoms.
(Top)
Obtaining a working initial PDB file:
Structure 3v03.pdb is a dimer. We select the second chain:
gmx pdb2gmx -f 3v03B_noH.pdb -o 3v03B.pdb -p 3v03B.top -i 3v03B.itp -water spce -ff oplsaa &> 3v03B_pdb2gmx.log
to obtain [3v03B.pdb](simulation/monomer/solvatedBSA/3v03B.pdb), among other output files.
Solvating the protein with water:
Use OPLSAA force field and SPCE water model:
gmx pdb2gmx -f 3v03B.pdb -o 3v03B.gro -p 3v03B.top -i 3v03B.itp -water spce -ff oplsaa &> 3v03B_pdb2gmx.log
outputs:
- 3v03B_pdb2gmx.log states 579 residues, 9173 atoms, and a total charge of -16e. We need to place 16 Na atoms in order to neutralize the system.
- 3v03B.gro processed binary gromacs file.
- 3v03B.itp positional restraints file.
- 3v03B.top topology of the isolated molecule.
The below command centers the protein in the box (-c), and places it at least 1.0 nm from the box edge (-d 1.0). The box type is defined as a cube (-bt cubic):
gmx editconf -f 3v03B.gro -o 3v03B_newbox.gro -c -d 1.0 -bt cubic
Fill with solvent (water). Solvation is accomplished using solvate:
cp 3v03B.top 3v03B_solv.top
gmx solvate -cp 3v03B_newbox.gro -cs spc216.gro -o 3v03B_solv.gro -p 3v03B_solv.top -scale 0.7
Add ions ([ions.mdp](simulation/monomer/solvatedBSA/ions.mdp)):
cp 3v03B_solv.top 3v03B_solv_ions.top
gmx grompp -f ions.mdp -c 3v03B_solv.gro -p 3v03B_solv_ions.top -o ions.tpr
gmx genion -s ions.tpr -o 3v03B_solv_ions.gro -p 3v03B_solv_ions.top -neutral # select group number 13.
Minimize the protein:
(In subdirectory simulation/monomer/solvatedBSA/minimize)
Minimize options specified in file [minimize.mdp](simulation/monomer/solvatedBSA/minimize/minimize.mdp).
gmx grompp -f minimize.mdp -c 3v03B_solv_ions.gro -p 3v03B_solv_ions.top -o em.tpr # creates binary input em.tpr.
gmx mdrun -ntmpi 16 -ntomp 1 -pin on -deffnm em
Relax the water around the protein:
(In subdirectory simulation/monomer/solvatedBSA/minimize/relax)
50ps NVT run using restraints on the protein atoms. Coordinates defining the equilibrium positions for the restraints are in file 3v03B.itp. This file is encoded in the topology file 3v03B_solv_ions.top. Options for the run specified in file [relax_1.mdp](simulation/monomer/solvatedBSA/minimize/relax/relax_1.mdp).
gmx grompp -f relax_1.mdp -c em.gro -p 3v03B_solv_ions.top -o relax_1.tpr
gmx mdrun -ntmpi 16 -ntomp 1 -pin on -deffnm relax_1 # use 16 MPI threads, pin each thread to a diffeent core.
Continue with a 500ps NPT run using restraints on the protein atoms. The goal is to shrink the simulation box to tipycal densities. Input file is [relax_2.mdp](simulation/monomer/solvatedBSA/minimize/relax/relax_2.mdp).
gmx grompp -f relax_2.mdp -c relax_1.gro -t relax_1.cpt -p 3v03B_solv_ions.top -o relax_2.tpr
time gmx mdrun -ntmpi 16 -ntomp 1 -pin on -deffnm relax_2
Recover some thermodynamics of the cell shrinkage:
gmx energy -f relax_2.edr -o relax_2.xvg
Center protein in the simulation box:
gmx trjconv -f relax_2.cpt -s relax_2.tpr -pbc mol -ur compact -center -o relax_2_centered.gro -boxcenter tric
Select group for centering: Group 1 ( Protein) has 9173 elements Select group for output: Group 0 ( System) has 84738 elements
Produces file relax_2_centered.gro.
Select closest N water molecules to protein:
The mass of the protein is 66054.497 a.m.u according to 3v03B_pdb2gmx.log. Thus a hydration level $h=0.1$ corresponds to 6605.4497 a.m.u of water, corresponding to 6605.4497/18.01528 = 366 water molecules.
In the experiment, the estimated hydration level is 0.17. Thus, we have to include 1.7*366 = 622 molecules.
Save the following selection to file selection_h0.17_waters.txt:
close_solvent_atoms = group "SOL" and within 0.1874 of group "Protein";
close_solvent_molecules = same resindex as close_solvent_atoms;
group "Protein" or close_solvent_molecules;
We found a cutoff = 0.1874 to retain the 622 closest water molecules, simply by counting the number of lines containing "SOL" in file stripped_system.pdb (see below).
Create index file stripped_system_h0.17.ndx with protein plus close waters, and also PDB file stripped_system_h0.17.pdb
gmx select -f relax_2_centered.gro -s relax_2_centered.gro -sf selection_h0.17_waters.txt -on stripped_system_h0.17.ndx -ofpdb stripped_system_h0.17.pdb -pdbatoms selected
(Top)
Silica $(SiO_2)_N$ has a sort-range order, whereby a Si atom is located at the center of a tetrahedron, whose vertices are defined by four O atoms.
In subdirectory silica/amorphous_from_md/ and also Materials Studio project silica in virtual machine LEIOAVIRTUAL.
I downloaded 21A box from a post in the Materials Studio community].
Inserting the protein in the silica box:
We start with file relax_2_centered.gro from directory simulation/monomer/solvatedBSA/minimize/relax, which is contained in a cubic box of 9.41445nm. We have to insert the protein in the silica box SiO2_107A_3d_NoConnect.pdb, which is a cubic box of 10.6974nm. Thus, we resize the box of the protein plus water:
gmx editconf -f relax_2_centered.gro -o relax_2_centered_resized.gro -bt cubic -box 10.6974
Next we select the 622 water molecules closest to the protein as we did before:
gmx select -f relax_2_centered_resized.gro -s relax_2_centered_resized.gro -sf selection_h0.17_waters.txt -on stripped_system_h0.17.ndx -ofpdb stripped_system_h0.17.pdb -pdbatoms selected
We load stripped_system_h0.17.pdb and SiO2_107A_3d_NoConnect.pdb in VMD to check the protein is in the middle of the silica box:
We need to "clean up" files stripped_system_h0.17.pdb and SiO2_107A_3d_NoConnect.pdb before merging protein and silica into a single PDB. We use merge_silica_and_protein.py for this purpose
python merge_silica_and_protein.py creates the following files:
Removing $SiO_2$ molecules overlapping the protein and water
(tbd)
(Top)
In subdirectory simulation/silica/cristobalite/cristobalite_2_2_2
Create topology with x2top utitily
We have modified oplsa.ff/ force field in simulation/silica/cristobalite/cristobalite_2_2_2/oplsa.ff/ accounting for silicon type SIO and oxygen type OSI present in the crystal. We modified files ffnonbonded.itp, ffbonded.ipt, and atomtypes.atp.
We create topology file with gmx x2top. This utitility requires description of connectivity between atom types in file oplsa.f/atomname2type.n2t.
gmx x2top -v -f cristobalite_2_2_2.gro -o cristobalite_2_2_2.top -ff oplsaa -name SIL -pbc -noparam
After this, a minimization (minin) followed by equilibration at either constant volume (relax_1) or constant pressure (relax_2) are carried out. Run relax_2 is continued with another NPT run (relax_3) which can be used to gather data. Script gmx_commands.sh containst the relevant gmx submissions.
Relaxation at constant pressure relax_2.mdp:
](simulation/silica/amorphous_from_md/confineBSA/relax_2_centered_resized.gif){kind=link}
](simulation/silica/amorphous_from_md/confineBSA/stripped_system_h0.17.gif){kind=link}
verlet-buffer-tolerance=-1 to use nstlist=1 and rlist. tau_p=100 for slow relaxation of the volumeContinuation at constant pressure relax_3.mdp:
gmx grompp -f relax_3.mdp -t relax_2.cpt -c relax_2.tpr -o relax_3.tpr # checkpoint requiredverlet-buffer-tolerance=-1 to use nstlist=10 (Top)
From the Crystallography Open Database we search for the crystal structure of alpha-cristobalite, one of the morphologies of $SiO_2$:
In Materials studio, the supercell replicating the unit cell 21 times along X and Y and 15 times along Z is able to surround BSA.
A post in materials studio prescribes a recipe to build amorphous silica starting with the amorphous cell module. We can take parts of this recipe to create amorphous silica from a cristobalite supercell.
From: Jason DeJoannis I think we built that structure in partnership with Corning. It is documented in this paper: Vyas, Dickinson and Armstrong-Poston, Mol Sim, 32, 2, 135 (2006). Build bulk glass with Amorphous Cell Pack cell with SiO2 and Si2O4 fragments at density of 2g/mL cvff_aug forcefield with Ewald Anneal the bulk NVT at 7000K for 20ps NPT for 20ps Reduce T in 1000K increments with 20ps at each temperature down to 300K Create the surface Cleave and add vacuum slab 50-60A Extend surface supercell as needed Anneal the surface Start at 600K and cool to 300K in steps of 100
(Top)
In Materials studio, the supercell replicating the unit cell 21 times along X and Y and 15 times along Z is able to surround BSA.
Inserting the protein in the silica box:
We start with file relax_2_centered.gro from directory simulation/monomer/solvatedBSA/minimize/relax, which is contained in a cubic box of 9.41445nm. We have to insert the protein in the cristobalite box cristobalite_21_21_15.pdb, which has dimensions 10.4406x10.4406x10.3835 nm. Thus, we resize the box of the protein plus water:
gmx editconf -f relax_2_centered.gro -o relax_2_centered_resized.gro -bt triclinic -box 10.4406 10.4406 10.3835
Next we select the 622 water molecules closest to the protein as we did before:
gmx select -f relax_2_centered_resized.gro -s relax_2_centered_resized.gro -sf selection_h0.17_waters.txt -on stripped_system_h0.17.ndx -ofpdb stripped_system_h0.17.pdb -pdbatoms selected
We need to "clean up" files stripped_system_h0.17.pdb and cristobalite_21_21_15.pdb before merging protein and silica into a single PDB. We use merge_silica_and_protein.py for this purpose.
python merge_silica_and_protein.py creates the following files:
Removing $SiO_2$ molecules overlapping the protein and water
The cristobalite lattice projected onto the XY plane can be viewed as a distorted squared lattice with the Si atoms occupying the corners, and the oxygen atoms located in the middle of the square sides. See below:
With this arrangement, we can univoquely associate two oxygens to each Si, out of the four first-neighbors. We chose the oxygen with the largest X-component of its position relative to the Si, and the oxygen with the largest Y-component. This way we can paritition the cristobalite crystal in a unique way.
Script carve_silica.py uses silicanet.py takes the previous partition to create a hole in the cristobalite lattice. It takes as input file SiO2_protein_0.17.pdb and partitions the cristobalite lattice into $SiO_2$ molecules. Molecules having any of its atoms overlapping the protein+water are removed. The overlapp parameter is tried with different values, from 1.0 to 2.5 Angstroms, producing different carved files. Also the names of the silica atoms will match those of the force field in subdirectory simulation/silica/cristobalite/confineBSA/poretop/confinedBSA.ff/:
From the PDB files, it seems that overlapp > 1.5 leaves too many $SiO_2$ molecules and steric clashes that will not be resolved by a minimization run.
Simulating the confining cristobalite without bonds
As a first approximation, we assume each cristobalite atom is independent of the others, and it's subjected to a harmonic restraint. For cristobalite atoms in the interface of the protein+water, we can use as spring constant the value extracted from the experimental MSD of the dry protein (~2E-06 nm^2/K). For the bulk cristobalite atoms located away from the protein+water, we can use the spring constant derived from the simulations of the small cristobalite (~7E-07 nm^2/K).
in directory simulation/silica/cristobalite/confineBSA/poretop/nobonds we have modified force field oplsaa.ff/ containing nonbonded interactions for the cristobalite atom types SIO, OSI, and OA. We modified the following files:
We use SiO2carved_ovl1.5_protein_0.17.pdb and create a suitable .gro file with script adaptPDB.py.
python adaptPDB.py # produces confinedBSA_0.gro
Create topology and restraint files:
gmx pdb2gmx -f confinedBSA_0.gro -o confinedBSA.gro -p confinedBSA.top -i confinedBSA.itp -water spce -ff oplsaa &> confinedBSA.log
produces the following files:
The total charge of the system is reported to be -342.900. To make the system neutral, we increase the positive charge of the SIO atoms. There are 23694 SIO atoms in the silica topology file confinedBSA_Other2.itp. Thus, the new charge will be 1.1 + 342.9/23694 = 1.1144720182324641. We edit file confinedBSA_Other2.itp and make this change.
We also updated the SIO charge in the following files of the force field:
We renamed Other2 to Silica, modifying confinedBSA.top and producing producing new files confinedBSA_Silica.itp and confinedBSA_Silica_pr.itp
Topology for the carved cristobalite
script evaluate_bonding.py uses silicanet.py to calculate average bond and angle distances in the cyrstobalite using the cristobalite in SiO2_protein_0.17.pdb. The results:
bond_Si_O = 1.57833041748
angle_O_Si_O = 109.149737561
angle_Si_O_Si = 146.498743845
Working directory: simulation/silica/cristobalite/confineBSA/poretop/yesbonds
We use Zhang's force field with the previously found bond and angle average quantities. Resulting force field in subdirectory simulation/silica/cristobalite/confineBSA/poretop/yesbonds/oplsaa.ff/.
From Simulating the confining cristobalite without bonds, we now that the charge of SIO should be 1.1144720182324641 to render a neutral silica+protein+water system. The charges are thus:
If we only want to simulate the silica, then a charge of SIO equal to 1.113796741791171 will render a neutral silica sytem.
To create a topology file for the carved cristobalite (without the protein) we edit SiO2carved_ovl1.5_protein_0.17.pdb, removing the protein and water atoms, and then saving as SiO2carved_ovl1.5.pdb. Then, script topology_silica.py will input this file and use module silicatop.py to generate an include topology file. we produce SiO2carved_ovl1.5_spw0.itp for a neutral silica+protein+water system, and SiO2carved_ovl1.5_s0.itp for a neutral silica system.
We simulate the silica alone in directory simulation/silica/cristobalite/confineBSA/poretop/yesbonds/onlysilica/ by copying SiO2carved_ovl1.5_s0.itp into SIL.itp. We had to tweak the last SIO with a slightly different charge (1.11505874179 instead of 1.11379674179) so that gmx grommp would report a neutral system. This is a numerical inaccuracy of gmx grommp. Also we had to comment the info in the [ defaults ] section because this info is previously set up in oplsaa.ff/forcefield.itp.
With script shift_indexes.py we shift the residue and atom indexes of SiO2carved_ovl1.5.pdb to create confiningSIL.pdb. This file starts with atom number one and has residue number equal to one.
We manually created confiningSIL.top and minimize.mdp
gmx editconf -f confiningSIL.pdb -o confiningSIL.gro
gmx grommp -f minimize.mdp -c confiningSIL.gro -p confiningSIL.top -o minimize.tpr
gmx mdrun -ntmpi 16 -ntomp 1 -pin on -deffnm minimize
The minimization movie shows an expanding hole, causing a deformation in the crystal structure close to the hole surface. The expansion is due to the net negative charge caused by the dangling oxygens. These oxygens repel each other, and they quickly shield from each other by "turning" towards the bulk crystal. There are some $SiO_2$ molecules in the hole flying around. The expansion due to the negative charge will be less pronounced once we include the protein because it will act as a shielding insulator, but the steric clashes between protein and surface with induce an expansion on the hole.
Inserting the protein in the carved crystobalite with bonds
In directory simulation/silica/cristobalite/confineBSA/poretop/yesbonds/wholesystem/ we copied SiO2carved_ovl1.5_spw0.itp into SIL.itp and modified this file as follows:
](simulation/silica/cristobalite/confineBSA/relax_2_centered_resized.gif){kind=link}
](simulation/silica/cristobalite/confineBSA/stripped_system_h0.17.gif){kind=link}
[ defaults ] sectiongmx grompp complained the system had no net zero charge.We copied file SiO2carved_ovl1.5_protein_0.17.pdb into confinedBSA.pdb
We created topology file confinedBSA.top by hand, retrieving all info for the protein and water from previous topology file of the solvated BSA, 3v03B_solv_ions.top
Input run files minimize.mdp and relax_1.mdp
gmx editconf -f confinedBSA.pdb -o confinedBSA.gro
gmx grompp -f minimize.mdp -c confinedBSA.gro -p confinedBSA.top -o minimize.tpr
gmx mdrun -ntmpi 16 -ntomp 1 -pin on -deffnm minimize
gmx grompp -f relax_1.mdp -c minimize.gro -p confinedBSA.top -o relax_1.tpr
gmx mdrun -ntmpi 16 -ntomp 1 -pin on -deffnm relax_1 # 4ns in the NVT
gmx grompp -f relax_2.mdp -c relax_1.gro -p confinedBSA.top -o relax_2.tpr
gmx mdrun -ntmpi 16 -ntomp 1 -pin on -deffnm relax_2 # 2ns in the NPT to relax the box size
gmx grompp -c relax_2.tpr -t relax_2.cpt -p confinedBSA.top -f equil_1.mdp -o equil_1.tpr # continuation but in the NVT
gmx mdrun -ntmpi 16 -ntomp 1 -pin on -deffnm equil_1
The 2ns relaxation at NPT (relax_2) proved insufficiently long since the pressure became negative in the equil_1 run. Thus, we extend the simulation for 10ns in Titan by appending 1ns jobs. For this we need to transfer to Titan all files relax_2.* and the following
ln -s relax_2.tpr relax_2.0.tprIn Titan, we just run the substitution script by passing as argument the number of the firs run:
./relax_2.submitnext.sh 0
Will trigger submission of the first nanosecond, followed by 9 more runs up to a total of 10ns
After relaxation, we run an equilibrium run for 1ns in the NVT ensemble, and then append subsequent runs for a total of 12ns:
ln -s equil_1.tpr equil_1.0.tpr and then start the chained submission scripts with ./equil_1.submitnext.sh 0 .(Top)
We simulate four systems in order to find out which system provides us with the best comparison to experimental MSD at 300K:
We already simulated the first system. The other three systems we create and simulate in subdirectories simulation/confinedBSA/h017_q0.55_ovl1.5, simulation/confinedBSA/h017_q0.55_ovl2.5, and simulation/confinedBSA/h017_q0.90_ovl2.5, respectively.
With gmx msd we calculate MSD on the equil_1.xtc trajectories
The system $q_{OA}=0.90$ and silica-system overlap 2.5A approaches the most the experimental point. We select this system for the temperature scan.
(Top)
Files confinedBSA.gro and confinedBSA.pdb are same as for system $q_{OA}=0.55$ and overlap=2.5A. We use topology_silica.py to create SiO2carved_q0.90_ovl2.5.itp which we rename to SIL.itp.
As usual, we ran gmx grompp -f minimize.mdp -c confinedBSA.gro -p confinedBSA.top -o minimize.tpr and use the reported excess charge to edit file SIL.itp and there tune the charge of the SIO atoms to render a neutral system.
We run the following jobs at 300K in Titan:
Temperatures scan within subdirectory simulation/confinedBSA/h017_q0.90_ovl2.5/Tscan
We run the temperatures scan in Titan from $T_i=40$ to $T_e=420$ in increments of $20K$. Preparation of input files prior to submitting jobs to the queue done by running bash_commands.sh.
The following jobs:
(Top)
https://bionanostructures.com/interface-md/
Latest version 1.5 in subdirectory simulation/forcefields/INTERFACE_FF_1_5/FORCE_FIELDS. Contains the force-field in PCFF and CHARMM format, plus a series of bulk and surface silica systems.
(Top)
Prof. Yingchun Liu sent me a force-field for silica after a personal request. This after reading Zheng16. This paper referenced Wensink00 paper.
Located in subdirectory simulation/forcefields/Zheng16. Contains bonded and nonbonded parameters.
(Top)
A possible series of steps:
(Top)
Quoted text
image caption |
some text |
In [ ]: