Difference between revisions of "How to do MDFF"
Line 62: | Line 62: | ||
==Analysis== | ==Analysis== | ||
===Check contacts of atoms=== | ===Check contacts of atoms=== | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
Basically it is using "measure contacts" function of VMD and list indexes of CA from the two selections and within 10 A. It is thus recommended to read only CA coordinates as your input. As we are interested in the residue ID of CA, the rest of the code then converts indexes into their corresponding resid. The resulting file will be two columns of integers listing pairs of CA which is within 10 A. | Basically it is using "measure contacts" function of VMD and list indexes of CA from the two selections and within 10 A. It is thus recommended to read only CA coordinates as your input. As we are interested in the residue ID of CA, the rest of the code then converts indexes into their corresponding resid. The resulting file will be two columns of integers listing pairs of CA which is within 10 A. | ||
In bash you may have a quick look at your result by | In bash you may have a quick look at your result by | ||
cat result.dat | cut -d ' ' -f1 | sort | uniq -c | cat result.dat | cut -d ' ' -f1 | sort | uniq -c | ||
− | in which you may change f1 to f2 as to view only column 2 of the file instead of only column 1 | + | in which you may change f1 to f2 as to view only column 2 of the file instead of only column 1. |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
===Scripts=== | ===Scripts=== |
Revision as of 23:39, 24 October 2014
Set-up
- Prepare initial structure
- !!!Be very careful to your initial structure!!!
- Autopsf
- MInimize your initial structure by first fixing backbone for 5000 and then no-fix for 10000 steps
- Solvate
- Add ion
- Prepare density map file (mrc to trimmed.dx)
- !!!Be very careful to the box sizes so that the water box must include the whole necessary density map!!!
- Prepare MDFF initial files
package require mdff
if {0} {
mdff gridpdb -psf ionized.psf -pdb ionized.pdb -o ionized-grid.pdb
package require ssrestraints
ssrestraints -psf ionized.psf -pdb ionized.pdb -o ionized-extrabonds.txt -hbonds
package require cispeptide
package require chirality
mol new ionized.psf
mol addfile ionized.pdb
cispeptide restrain -o ionized-extrabonds-cispeptide.txt
chirality restrain -o ionized-extrabonds-chirality.txt
}
- Prepare MDFF namd files - the first is for heating (manually change ITEMP to 0 and FTEMP to 310)
mdff setup -pbc -o ./output/symmetry-2-heat -psf ionized.psf -pdb ionized.pdb -griddx grid-trimmed.dx -gridpdb ionized-grid.pdb -extrab {ionized-extrabonds.txt ionized-extrabonds-cispeptide.txt ionized-extrabonds-chirality.txt} -gscale 0.1 -minsteps 10000 -numsteps 100000
for {set i 1} {$i < 10} {incr i} {
mdff setup -pbc -o ./output/symmetry-2-long -psf ionized.psf -pdb ionized.pdb -griddx grid-trimmed.dx -gridpdb ionized-grid.pdb -extrab {ionized-extrabonds.txt ionized-extrabonds-cispeptide.txt ionized-extrabonds-chirality.txt} -gscale 0.1 -numsteps 1000000 -step $i
}
- Create 2-symmetry.pdb and 2-matrices.txt
- Add symmetrical restraints to heat-step1 (manually)
- Add symmetrical restraints to parameters (by script)
for (( i=1; i<=1; i++ ))
do
sed -i '48,$d' symmetry-2-long-step$i.namd
cat symmetry.txt >> symmetry-2-long-step$i.namd
echo -e "\nsource mdff_template.namd" >> symmetry-2-long-step$i.namd
done
- symmetry.txt are as follows:
symmetryRestraints on
symmetryfile 2-symmetry.pdb
symmetryk 200
symmetryMatrixFile 2-matrices.txt
symmetryfirststep 201
symmetryfirstfullstep 1000200
- !!!Manually change ITEMP to 0 and FTEMP to 310 in heat-step1!!!
- !!!Manually change both TEMP to 310 in long-steps!!!
sed -i 's/TEMP 300/TEMP 310/g' symmetry-2-long-step*.namd
- gosub
Analysis
Check contacts of atoms
Basically it is using "measure contacts" function of VMD and list indexes of CA from the two selections and within 10 A. It is thus recommended to read only CA coordinates as your input. As we are interested in the residue ID of CA, the rest of the code then converts indexes into their corresponding resid. The resulting file will be two columns of integers listing pairs of CA which is within 10 A.
In bash you may have a quick look at your result by
cat result.dat | cut -d ' ' -f1 | sort | uniq -c
in which you may change f1 to f2 as to view only column 2 of the file instead of only column 1.
Scripts
Scripts cron_ana.sh has been written to run three tcl analysis scripts writepdb.tcl, cal_contacts.tcl and cal_dcd.tcl and some more modifications.
MATLAB
Now you have two columns of residue IDs and be ready to turn them into a histogram.
In addition to the histplot function in MATLAB, you may also use histc (hist-count) to do the sorting and counting first:
c = histc(VarName, (min, max));
for example:
c = histc(bot_sym_heat78,(4:255));
Then plot the variable c which ranges in 4:255.