Difference between revisions of "How to do MDFF"

Revision as of 22:25, 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

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. cal_contacts.tcl

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.