Difference between revisions of "How to do MDFF"
m |
m |
||
| Line 1: | Line 1: | ||
| + | ==Set-up== | ||
| + | *Prepare initial structure | ||
| + | *'''!!!Be very careful to your initial structure!!!''' | ||
| + | *Autopsf | ||
| + | *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 | ||
| + | <code> | ||
| + | package require mdff<br /> | ||
| + | if {0} {<br /> | ||
| + | mdff gridpdb -psf ionized.psf -pdb ionized.pdb -o ionized-grid.pdb | ||
| + | |||
| + | package require ssrestraints<br /> | ||
| + | ssrestraints -psf ionized.psf -pdb ionized.pdb -o ionized-extrabonds.txt -hbonds | ||
| + | |||
| + | package require cispeptide<br /> | ||
| + | package require chirality<br /> | ||
| + | mol new ionized.psf<br /> | ||
| + | mol addfile ionized.pdb<br /> | ||
| + | cispeptide restrain -o ionized-extrabonds-cispeptide.txt<br /> | ||
| + | chirality restrain -o ionized-extrabonds-chirality.txt<br /> | ||
| + | } | ||
| + | </code> | ||
| + | *Prepare MDFF namd files - the first is for heating (manually change ITEMP to 0 and FTEMP to 310) | ||
| + | <code> | ||
| + | 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} {<br /> | ||
| + | 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<br /> | ||
| + | } | ||
| + | </code> | ||
| + | *'''!!!Manually change TEMP to 310!!!''' | ||
| + | <code> | ||
| + | sed -i 's/TEMP 300/TEMP 310/g' symmetry-2-long-step*.namd | ||
| + | </code> | ||
| + | *'''!!!Manually change TEMP to 310!!!''' | ||
| + | <code> | ||
| + | </code> | ||
| + | |||
==Analysis== | ==Analysis== | ||
===Check contacts of atoms=== | ===Check contacts of atoms=== | ||
Revision as of 15:51, 22 October 2014
Contents
[hide]Set-up
- Prepare initial structure
- !!!Be very careful to your initial structure!!!
- Autopsf
- 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
}
- !!!Manually change TEMP to 310!!!
sed -i 's/TEMP 300/TEMP 310/g' symmetry-2-long-step*.namd
- !!!Manually change TEMP to 310!!!
Analysis
Check contacts of atoms
mol new solvent_CA_last.pdb
set name1 "P3 P4"
set name2 "P7 P8"
set sel1 [atomselect top "segname $name1"]
set sel2 [atomselect top "segname $name2"]
set res [measure contacts 10 $sel1 $sel2]
set col1 [lindex $res 0]
set col2 [lindex $res 1]
set outf [open $name1-$name2.txt "w"]
for {set i 0} {$i < [llength $col1]} {incr i} {
set num1 [lindex $col1 $i]
set num2 [lindex $col2 $i]
set sel3 [atomselect top "index $num1"]
set out1 [$sel3 get resid]
set out3 [$sel3 get segname]
set sel4 [atomselect top "index $num2"]
set out2 [$sel4 get resid]
set out4 [$sel4 get segname]
puts "$out1 $out2"
puts $outf "$out1 $out2 $out3 $out4"
}
quit
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.
A check-XXX.dat is also written with segname of the pairs of CA as to check whether they come from the same segname or not. Also use
cat result.dat | cut -d ' ' -f1 | sort | uniq -c
to check their repetition.
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.
