Setup and equilibration of protein

Here focus is on the general setup (adding hydrogens and determining protonation state) with the maestro program. We additionally go through how maestro can be checked by Ulf Rydes programs. A complete list of the setup can be found at Ulf Ryde's page. The procedure below goes through steps 1-7.

Table of Contents

1. Characterization of the pdb
2. Protein setup
3. Non-standard residues
4. Rename residues
5. Add non polar hydrogens, missing atoms and solvent cap with tleap
6. Structure check and optional charge modification
7. Force field check by energy calculation
8. Simulated annealing
9. Write pdb file

1. Characterization of the pdb

Check (use open-office chart):

Is the protein a monomer, dimer ? or n-mer?
Are there Cys-Cys-crosslinks (see also 4. Rename residues below).
Are there missing parts of the structure (the first amino acids are often not resolved, check REMARK 465)? We may later add missing atoms automatically using tleap.
Are the non-standard residues? Note down the residue numbers and names. Check if there a force field available for these residues.
At what pH was the structure obtained?

2. Protein setup

Protein setup with Maestro

Make log-file (write down important information concerning flipped residues etc.). Maestro can write a command log file automatically: “Edit” => “Preferences...” => “General”, submenu “Quitting Maestro” => Check “Write command log file”.
Open Prep. Wiz. and "Import" the pdb file.
Make sure "Delete water beyond 5 Å" is unchecked and hit "preprocess".
Review "Problems => Current overlapping atoms" If all close contacts are hydrogens, these will be optimized later and are likely not a problem in the end. In case they are water molecules, we return to them in 7. Force field check by energy calculation).
Review the "Problems => ".
- “alternative positions". By "Next position" you can see the other position(s). By clicking "commit" you select one configuration, keep those with the highest occupancy (if close to active side, consider to calculate both). Mark down which you choose and why (with 50:50 it is basically guessing).
- “overlapping atoms". If caused by water molecules, check again after the Simulated annealing (step 8.) below. If the problems are not resolved by the optimization, then note down the problematic residues/atoms and get back to it in step 7.
- “missing atoms” If parts of residues are missing, note the residues down and check again visually after tleap filled the missing parts (see Step 5).
Go to "Review and Modify" and then "Analyze Workspace" and delete undesired parts (e.g. co-crystallized solvent molecules that are not water or other artefacts).
Save a pdb before optimization.
Optimize H-positions under "Refine" (can take some minutes; usually 5-10 but might take up to several hours for large systems).
Run Propka and check his residues (see also below).
- Check His (normally close to 7). We will deal with His in more detail below.
- Check if Asp and Glu are above 7 (they are then protonated, otherwise not).
- Chek of Lys, Arg, and Tyr have pKa below 7 (then they are not protonated).
- Note propka is a good tool, but not always reliable. Use it as guidelines (pKa calculations are very difficult).
Check histidines with changepdb (see Analysis with changepdb below) - use the pdb from the pdb database.
Check burried charges through visualization and with changepdb (see Analysis with changepdb below). There should be none that are not involved in ionic pairs - excpetions for metal ligands.
Decide on protonation state for Asp, Glu, Lys, Arg, Tyr and His
Use "Interactive optimizer => Analyze network"
Note (and check) which Asn and Gln groups that are flipped. They can be flipped by hitting the arrow to the left/right side in Maestro.
Check His protonation states (HIE, HID or HIP). Does it fit with analysis from changepdb? If not, then change them by hitting the arrow to the left side.
Check that all residues that should be charged are actually charged.
Lock the histidines by checking the "lock" box. Protonation states as well as the state of other (e.g. flipped) residues that should remain as they are should also be checked as "locked". Hit "Optimize".
For metal enzymes, it can be a good idea to check the metal site after the optimization.
Save a pdbfile (HID, HIE, HIP must be specified afterwards - this can be done with changepdb command "hi").

Analysis with changepdb

Histidine analysis
- run changepdb.
- give the input pdb.
- command: hh (creates output to screen and files "hisX" with X = 1,2,3...).
- copy the output to screen into a file (e.g. named "hishb").
- grep after the relevant information (grep HHb hishb).
- inspect the histidine files and decide on protonation state (HIE, HID, HIP). Inspect the hydrogen-bonding network around the histidines - for highly negatively charged proteins, solvent exposed histidines are generally HIP (unless there are good arguments against it). Their pKa values can also be calculated with PROPKA. If the histidine is surrounded by many negatively charged residues, it is more likely that it is protonated. Also look at the hydrogen bond network to decide whether HIE or HID are more likely (keep in mind that some residues can be hydrogen bond donors and acceptors). Overall HIE is more common than HID.
Charge of the complete system
- Run changepdb
- give the input pdb
- command cc
- Use charges
  - "Y" requires a pdb file with charges (they are not in crystal structures)
  - "N" can run with the crystal structure pdb.
Burried charges
- run changepdb.
- give the input pdb.
- command bc.
Cys-Cys cross-links
- run changepdb
- give the input pdb
- command cy
Metal centers (typically not done)
- run changepdb.
- give the input pdb.
- command: des (gives file *.actX (X = 1,2,..) that contains an active site cut-out.
Change histidines
- run changepdb.
- give the input pdb.
- run command "hi".
- type "D" for HID, "E" for HIE, "P" for HIP.
- type "w" (save) and then "q" (quit).

3. Non-standard residues

Non-standard residues requires determination of charges and sometimes a "dummy" force field. More advanced force fields can be constructed if necessary (see PON parametrization)

Cut-out the non-standard residue. If it is covalently bound as many metal active sites, see Cut out a covalently bound residue
Run a geometry optimization on the non-standard residues [typically B3LYP/def2-SV(P) or TPSS/def2-SV(P)]. Make sure that the structure is not too different from the crystal structure or other trustworthy experimental structures (if there are large deviations, consider using a single point calculation or only optimized hydrogens to generate charges).
- When only optimizing H-atoms: add f behind all non-hydrogen atoms in coord and run define (use change coordinates? –> yes –> ired).
Insert "$esp_fit resp" in the turbomole "control" file and run turbomole (must be version > 7.0):
```
ridft -resp -proper > logm
```
The use of "-proper" indicates that no wave function optimization is required. Delete it if you did not do a wave function optimization. Look for "GRID" in logm to check whether the Electrostatic Potential (which will be used to fit the resp charge) is generated.

Note: It can be convenient to ensure that the order of atoms is the same as in the pdb file. This can be ensured by using short-cut to generate "control" and "coord" files for turbomole. Make a syst1 file and use pdbtocomqum followed by comqumtoturbo as described under "Setup of ComQum Calculations" in the ComQum section.

Use now the script "changepot" (generates resp.pot, resp.in1, resp.in and qin). Note: changepot might give an error due to a missing line in logm. Insert this line with:
```
awk '/charge for group # 1/{x++} x==2{sub(/charge for group # 1/,"writing ESP data for Amber to file esp.dat \n&")}1' logm > logm_2 && mv logm_2 logm
```
changepot:
1. Select program, type t for turbomole.
2. Press enter (new); this is the name of the potential file (new is default).
3. Name of turbomole output (logm).
4. Do not use Boltzman weights (press enter).
5. Press enter (special requirements if # grid pints > 99999).
6. Press enter (output: resp.pot).
7. Press enter (new format)
8. Press enter (resp.in / resp.in1).
9. Press enter (qin).
  
  The files resp.in and resp.in1 are input files for the fitting of resp charges. qin are initial charges. Fitting of resp charges is usually done in two steps:
  1. First a fit involving all atoms with loose convergence. This is done with resp.in.
  2. In a next step, all heavy atoms are fixed while hydrogens in methylene (CH2) and methyl (CH3) groups are constrained to be identical and the fit is done with more tight convergence.
  Steps (i) and (ii) are done in the next few points below:

Check if there are any polar hydrogens bound to the same atom (e.g. water). These should be constrained to have the same charge (this is done in resp.in, see below for the format):

Format of resp.in files:

&cntrl  ihfree=0, qwt=0.0005, iqopt=2 # qwt = optimization parameter
&end
1.0
Mol1
2  63 # charge and tot. # of atoms
7   0 # Atom nr. (here 7 = N) and command (-1 = freeze, 0 = optimized freely, n = fix to the same as atom #n e.g. 43 means fix to same value as atom number 43)
6    0
6    0

Run resp (part of the AMBER suite)

resp -O -i resp.in -o resp.out -q qin -e resp.pot

Open resp.in1: Fix all charges except carbon atoms with more than one hydrogen (usually CH3 or CH2 groups). Fix the the hydrogens to have the same charge. qwt should be 0.001. Save qout in qin1:
```
cp qout qin1
```
and run
```
resp -O -i resp.in1 -o resp.out1 -q qin1 -e resp.pot
```
The final charges are now in resp.out1 - add this charges to the pdb file you want to start the equillibrium from (see more below and under ### 6. Step 6.).
- Guide to adding charges to pdb: Column q(init) in resp.out1 gives the input charges. Column q(opt) are the output charges, which are the final charges. For the QM atoms (but not for the junction atoms), they replace the charges in the last column of the pdb file (make sure to copy the charges in the right order to the atoms). Note that with this replacement, the total charge of the protein will no longer be integer (corrected with junctions below).
- Junction atoms (atoms converted to a hydrogen atom in the quantum system): After adding charges to the QM atoms, it is likely that the total charge of the protein no longer correctly adds up (can be checked with changepdb, see next point). This is fixed by adding a uniform constant to all junction atoms (add this factor to the existing AMBER charge), ensuring that the charge adds up correctly.
Check the charge with changepdb (cc) - see Analysis with changepdb. Ensure that the charges add up.

Cut-out of covalently-bound residue

To extract pdb file and cap with hydrogens there are two options:

Use Marlisa’s python script to generate a capped xyz file. The xyz file can be transfered to tmole coord file format with openbabel (copy the pdb atom names into the coord file and note down the junction atoms).
Use the input-pdb from Step 5. below. For this pdb, use the scripts “pdbtocomqum” and “comqumtoturbo” to generate the turbomole coord and controle files. Note: “pdbtocomqum” and “comqumtoturbo” require a syst1 file. This syst1 file should only contain the residues you need RESP charges for. The format of the syst1 file can be found under "Setup of ComQum Calculations" in the ComQum section.

4. Rename residues

Cys => Cym if Cys is a metal ligand.
Cys => CYX if Cys is a Cys-Cys cross link.
HIS => HID, HIE or HIP if Maestro was used for the setup.
HOH => WAM if water is a metal ligand. Move these waters to the top of the list of water molecules and add TER before and after them. Note: that wam.in and wam.dat will be needed for tleap (see next point, below).
If pdb still contains ANISOU or CONNECT keyword, just remove it, for example with grep -v ANISOU inp.pdb > out.pdb

5. Add non polar hydrogens, missing atoms and solvent cap with tleap

Note: We here run the equilibration with an octahedral solvent cap, whereas older procedures used a spherical cap (less efficient). The old procedure and corresponding AMBER inputs can be found here.

Solvent cap

Data required for solvating the protein can be obtained through changepdb:

changepdb
1. run changepdb
2. Press c (cap)
3. Press m (move)
4. Press w (write)

Run tleap (see below for sample leap.in file)

tleap -f leap.in

Sample tleap file

source leaprc.ff14SB # Amber protein force field
loadAmberPrep hic.in # special residue prep.in file
loadAmberPrep cu2.in # special residue prep.in file
loadAmberParams hic.dat # special residue data file
loadAmberParams cu2.dat # special residue data file
x=loadpdb lpmo_setup_part_3_NH2corrected.pdb # pdb input file
bond x.56.SG x.178.SG # Cys-Cys cross-link
solvateOct x TIP3PBOX 10 # solvation (see above)
saveamberparm x prmtop prmcrd
quit

Note: Residue numbers for the Cys-Cys cross-link refer to the residue numbers in the pdb input and need to be updated

The newly generated prmtop and prmcrd files will be used for the equilibration. prmtop contains anything that does not change during the simulation, prmcrd contains the coordinates that will be changed during the simulation. The coordinates in prmcrd will be translated compared to the ones in the input pdb file.

There will likely be errors from tleap. A few examples on how to fix these are given below:

Typical errors in tleap

End-groups (OXT and AMT): If one wishes to have the end amino group not charged, then delete the lines in addPdbResMap. Copy file leaprc.ff03 to the local directory;
```
cp $AMBERHOME/dat/leap/cmd/leaprc.ff03 .
```
Alternatively, clearPdbResMap can be added to the tleap input file.

Prep.in files: For special residues, a prep.in file must be constructed for the atom-types. A pdb file with the residue can be cut out from the pdbfile. The prep file can be constructed by (amber)

antechamber -i hic.pdb -fi pdb -o hic.in -fo prepi -nc 0 -c bcc -rn HIC -at amber # AMBER force field
antechamber -i hic.pdb -fi pdb -o hic.in -fo prepi -nc 0 -c bcc -rn HIC -at gaff # GAFF force field

Missing parameters: In some cases, neither the regular amber protein force field nor the GAFF force field have parameters. In this case, the parameters must be constructed.
Missing atoms: tleap can usually add missing atoms, but be aware that it might be required to delete hydrogens that maestro added (to ensure tleap recognize the residue correctly). Check the error/warning from tleap fo see of this is the case
Missing residues: Many pdbs have missing residues, usually since they were too difficult to resolve crystallographically (often seen in flexible regions). Usually we do not fill them in, but Maestro can do this.
- If the last 10-20 residues are missing , we usually do not fill them in.
- If the residues are in the middle of the protein, they might need to be filled (important for MD simulations). If not, the issues in tleap can be resolved by adding TER after the first unresolved residue. To get a chargedlast residue but no charged groups for the missing residues use the following protocol:
  1. copy leaprc.protein.ff14SB or leaprc.ff03 to the working directoy as above (see under End-groups). Note: The altered leaprc.protein.ff14SB or leaprc.ff03 file needs to always be copied when starting QM/MM calculations (or its full path needs to be given in leap.in).
  2. delete everything in leaprc.protein.ff14SB or leaprc.ff03, except for the residue type in the end that should be charged.
Unperturbed charge is not 0 warning. This is related to Ewald summation, which needs a neutral system. Can be ignored (if possible, adjust surface histidine residues to get a zero net chage).

6. Structure check and optional charge modification

Build a pdb file by typing mkpdbout and rename the output (pdbout) to a new file (e.g. pdbout mv start.pdb). Inspect start.pdb (is the solvent cap reasonably placed? Is the structure looks reasonable?).
Check that the charges add up to an integer by using changepdb, option CC (see Analysis with changepdb. Note: For QM residues/atoms the RESP charges from should be in the start.pdb file. For the junctions, we use a modified charge (see last part of Step 3.).
After ensuring that the charges in start.pdb are the same as the ones added last in Step 3., the prmtop needs to be updated with these charges. Use the program changeparm:

changeparm <<EOF
prmtop
m
start.pdb
w

q
EOF

7. Force field check by energy calculation

It is a good idea to check the resulting energies with the program changeparm

changeparm
1. run changeparm
2. Select prmtop
3. Command e (for energy)
4. Command m (for mdrest)
5. Select prmcrd
- Output (first): 10 worst energies for bond, angle and dihedrals => Energies above 100 kcal/mol should be checked.
- Output (second): Lennard-Jones and electrostatic energies => We typically ignore the electrostatic, but vdW energies above 100 kcal/mol should be checked (over 1000 kcal/mol something is wrong):
  - If it is a crystal water, it can be deleted before Step 5. Repeat the step with this water deleted
  - Can also be due to the non-hydrogen atoms added by tleap. Can be solved with an MD run including only these residues.

8. Simulated annealing

To equilibrate the system, we use a procedure of 5 steps. The input files and a script to run all the calculations are given below.

mpirun pmemd.MPI  -O -i   sander.in00 -o sander.out00 -p prmtop -c prmcrd  -r mdrest00 -ref prmcrd # step 1
mpirun pmemd.MPI  -O -i   sander.in0 -o sander.out0 -p prmtop -c mdrest00 -r mdrest0 -ref mdrest00 # step 2
mpirun pmemd.MPI  -O -i   sander.in1 -o sander.out1 -p prmtop -c mdrest0 -r mdrest1 -ref mdrest0   # step 3
mpirun pmemd.MPI  -O -i   sander.in2 -o sander.out2 -p prmtop -c mdrest1 -r mdrest2 -ref mdrest1   # step 4
mpirun pmemd.MPI  -O -i   sander.in3 -o sander.out3 -p prmtop -c mdrest2 -r mdrest3 -ref mdrest2   # step 5

sander input files are given below. Note: All 1111 should be replaced with the last crystal water.
In the sander files we use restraint_wt=10000 (typical gives average movement of atoms by 0.003 A and 0.006 A maximum). Deviations in metal-ligand distances of up to 0.01 A. Setting restraint_wt=1000 typically gives too large average movement: 0.03 A on average and 0.3 A maximum. Deviations in metal-ligand distances of up to 0.13 A.

sander.in00 (Minimization without shake).

Minimization without shake
&cntrl
irest=0,ntx=1,
nstlim=0,dt=0.0005,
imin=1,maxcyc=1000,drms=0.001,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=1,ntf=1,
nsnb=40,cut=8.0,dielc=1.0,
ntpr=100,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
ntb=1,ntp=0,taup=0.2,
ipol=0,igb=0,
ncyc=10,ntmin=1,dx0=0.01,
ntr=1,restraint_wt=10000,
restraintmask=":1-1111 & !@H=" # replace 1111 with last crystal water
&end

sander.in0 (10 ps MD simulation without shake)

10 ps MD simulation without shake
&cntrl
irest=0,ntx=1,
nstlim=20000,dt=0.0005,
imin=0,maxcyc=100,drms=0.001,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=1,ntf=1,
nsnb=40,cut=8.0,dielc=1.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
ntb=1,ntp=0,taup=0.2,
ipol=0,igb=0,
ncyc=10,ntmin=1,dx0=0.01,
ntr=1,restraint_wt=10000, 
restraintmask=":1-1111 & !@H=" # replace 1111 with last crystal water 
&end

sander.in1 (1 ns equilibration with shake and constant volume; NVT)

&cntrl
irest=1,ntx=5,
nstlim=500000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=2,
nsnb=15,cut=8.0,dielc=1.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
ntb=1,ntp=0,taup=0.2,
ipol=0,igb=0,
ntr=1,restraint_wt=10000,
restraintmask=":1-1111 & !@H=" # replace 1111 with last crystal water  
&end

sander.in2 (10 ns simulated annealing at constant pressure)

10 ns simulated annealing at constant pressure
&cntrl
irest=1,ntx=5,
nstlim=5000000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=2,
nsnb=15,cut=8.0,dielc=1.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
nmropt=1,
ntr=1,restraint_wt=10000,
restraintmask=":1-1111 & !@H=" # replace 1111 with last crystal water 
&end

&wt type='TEMP0',istep1=     0,istep2=1600000,value1=370.00,value2=370.00  &end

&wt type='TEMP0',istep1=1600001,istep2=5000000,value1=370.00,value2=  0.00  &end

&wt type='TAUTP',istep1=     0,istep2=2000000,value1=  0.20,value2=  0.20  &end

&wt type='TAUTP',istep1=200001,istep2=3200000,value1=  1.00,value2=  1.00  &end

&wt type='TAUTP',istep1=3200001,istep2=4000000,value1=  0.50,value2=  0.50  &end

&wt type='TAUTP',istep1=4000001,istep2=5000000,value1=  0.05,value2=  0.05  &end

&wt type='END' &end

&rst iat=0 &end

sander.in3: Minimisation with shake

Minimisation with shake
&cntrl
irest=0,ntx=1,
nstlim=0,dt=0.002,
imin=1,maxcyc=10000,drms=0.001,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=2,
nsnb=15,cut=8.0,dielc=1.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ncyc=10,ntmin=1,dx0=0.01,
ntr=1,restraint_wt=10000,
restraintmask=":1-1111 & !@H=" # replace 1111 with last crystal water 
&end

Known problems

Gas bubbles. Use a larger box of water molecules (larger value for solvateOct).
In the case that SHAKE complains, it can often be relieved by running a short minimization followed and short MD without SHAKE (see input "sander.in00" and sander.in0) and then continue the annealing with SHAKE. If the problem persists, then continue without SHAKE (note that the inputs below are for the CPU version and thus will maybe work better with the old procedure with spherical solvent cap - see Step. 5 under old methods here).
```
mpirun -bind-to core -np 20 sander.MPI -O -i sander.in1_no_shake -o sander.out1_no_shake -r mdrest1_no_shake -e mden1_no_shake -p prmtop -c prmcrd
mpirun -bind-to core -np 20 sander.MPI -O -i sander.in2_no_shake -o sander.out2_no_shake -r mdrest2_no_shake -e mden2_no_shake -p prmtop -c mdrest1_no_shake
```
The sander.in1_no_shake and sander.in2_no_shake are given below

Simulated annealing without shake (sander.in1_no_shake)

Title
 &cntrl
  irest=0,ntx=1,
  nstlim=1200000,dt=0.0005,
  temp0=300.0,ntt=1,tautp=0.2,
  ntc=1,ntf=1,
  nsnb=15,cut=15.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
  ntb=0,ntp=0,taup=0.2,
  ipol=0,igb=0,
  vlimit=20.0,nmropt=1,
  ibelly=1,bellymask=":WAT | @H="
 &end

 &wt type='TEMP0',istep1=     0,istep2= 400000,value1=370.00,value2=370.00  &end

 &wt type='TEMP0',istep1= 400001,istep2=1200000,value1=370.00,value2=  0.00  &end

 &wt type='TAUTP',istep1=     0,istep2= 400000,value1=  0.20,value2=  0.20  &end

 &wt type='TAUTP',istep1= 400001,istep2= 800000,value1=  1.00,value2=  1.00  &end

 &wt type='TAUTP',istep1= 800001,istep2=1000000,value1=  0.50,value2=  0.50  &end

 &wt type='TAUTP',istep1=1000001,istep2=1200000,value1=  0.05,value2=  0.05  &end

 &wt type='END' &end

 &rst iat=0 &end

Minimization (CHECK THIS) (sander.in2_no_shake)

Title
 &cntrl
  irest=0,ntx=1,
  nstlim=0,dt=0.0005,
  imin=1,maxcyc=10000,drms=0.001,
  temp0=300.0,ntt=1,tautp=0.2,
  ntc=1,ntf=1,
  nsnb=15,cut=15.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
  ntb=0,ntp=0,taup=0.2
  ipol=0,igb=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ibelly=1,bellymask=":WAT | @H="
 &end

9. Write pdb file

When all the sander jobs are completed, build an initial pdb file from the final coordinates using changeparm
```
changeparm <<EOF
prmtop
p
pdb3_equil.pdb
m
mdrest3
q
EOF
```
Characterise the pdb with changepdb:
1. Number of atoms and residues
2. Total charge, now including the actual charges (command cc)
3. The CAp centre (the centre of the protein, including waters) and the radius (in output under max dist). Command CA
4. Get the central residue of the protein, using
  1. command RW
  2. radius of 0 Ångstrom
  3. command CAP

run the program cpptraj with cpptraj prmtop ptraj.in. Use the ptraj.in file below ptraj.in

trajin mdrest3
trajout mdrest3-wrap restart
image origin center familiar com :98 # central residue found with changepdb
go

Run changeparm
```
changeparm <<EOF
prmtop
p
pdb3.pdb
m
mdrest3-wrap
q
EOF
```
pdb3.pdb is the final pdb file with an octahedral cap.

Depreciated protocol with spherical solvent cap

5. Add non polar hydrogens, missing atoms and solvent cap with tleap

For a spherical cap, we run tleap with a slightly different input file, given below:

Run tleap (see below for sample leap.in file)

tleap -f leap.in

Sample tleap file

source leaprc.ff14SB # Amber protein force field
loadAmberPrep hic.in # special residue prep.in file
loadAmberPrep cu2.in # special residue prep.in file
loadAmberParams hic.dat # special residue data file
loadAmberParams cu2.dat # special residue data file
x=loadpdb lpmo_setup_part_3_NH2corrected.pdb # pdb input file
bond x.56.SG x.178.SG # Cys-Cys cross-link
solvatecap x TIP3PBOX {0.0 0.0 0.0} 40 # solvation (see above)
saveamberparm x prmtop prmcrd
quit

8. Simulated annealing

Run a simulated annealing run with shake, using the input:

Title
&cntrl
irest=0,ntx=1,
nstlim=1200000,dt=0.002,
temp0=300.0,ntt=1,tautp=0.2,
ntc=2,ntf=2,
nsnb=15,cut=15.0,dielc=1.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
ntb=0,ntp=0,taup=0.2,
ipol=0,igb=0,
vlimit=20.0,nmropt=1,
ibelly=1,bellymask=":WAT | @H="
&end

&wt type='TEMP0',istep1=     0,istep2= 400000,value1=370.00,value2=370.00  &end
&wt type='TEMP0',istep1= 400001,istep2=1200000,value1=370.00,value2=  0.00  &end
&wt type='TAUTP',istep1=     0,istep2= 400000,value1=  0.20,value2=  0.20  &end
&wt type='TAUTP',istep1= 400001,istep2= 800000,value1=  1.00,value2=  1.00  &end
&wt type='TAUTP',istep1= 800001,istep2=1000000,value1=  0.50,value2=  0.50  &end
&wt type='TAUTP',istep1=1000001,istep2=1200000,value1=  0.05,value2=  0.05  &end

&wt type='END' &end

&rst iat=0 &end

After this, run a minimization with 10000 steps

Title
&cntrl
irest=0,ntx=1,
nstlim=0,dt=0.002,
imin=1,maxcyc=10000,drms=0.001,
temp0=300.0,ntt=1,tautp=0.2, <tt>  ntc=2,ntf=2,
nsnb=15,cut=15.0,dielc=1.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
ntb=0,ntp=0,taup=0.2
ipol=0,igb=0,
ncyc=10,ntmin=1,dx0=0.01,
ibelly=1,bellymask=":WAT | @H="
&end

Protein setup

Setup and equilibration of protein

1. Characterization of the pdb

2. Protein setup

Protein setup with Maestro

Analysis with changepdb

3. Non-standard residues

Cut-out of covalently-bound residue

4. Rename residues

5. Add non polar hydrogens, missing atoms and solvent cap with tleap

Typical errors in tleap

6. Structure check and optional charge modification

7. Force field check by energy calculation

8. Simulated annealing

Known problems

9. Write pdb file

Depreciated protocol with spherical solvent cap

5. Add non polar hydrogens, missing atoms and solvent cap with tleap

8. Simulated annealing

results matching ""

No results matching ""