24小时热门版块排行榜

返回列表

nufang19a

金虫 (正式写手)

应助: 0 (幼儿园)
金币: 833.9
散金: 103
红花: 1
帖子: 449
在线: 101.7小时
虫号: 968803
注册: 2010-03-12
性别: GG
专业: 分子模拟

[交流] 【讨论】大家看看我的RMSD图，给些指点，谢谢了已有3人参与

跑完平衡动力学，在求RMSD for individual residues时发现输出的RMSD图有些诡异，我用的是这个程序：
# % $Id: residue_rmsd.tcl,v 1.4 2006/03/06 23:56:46 timisgro Exp $
proc rmsd_residue_over_time {{mol top} res} {

# use frame 0 for the reference
set reference [atomselect $mol "protein" frame 0]
# the frame being compared
set compare [atomselect $mol "protein"]
#make a selection with all atoms
set all [atomselect top all]
#get the number of frames
set num_steps [molinfo $mol get numframes]
#open file for writing
set fil [open residue_rmsd.dat w]

foreach r $res {
set rmsd($r) 0
}

#loop over all frames in the trajectory
for {set frame 0} {$frame < $num_steps} {incr frame} {
puts "Calculating rmsd for frame $frame ..."
# get the correct frame
$compare frame $frame
      $all frame $frame
# compute the transformation
set trans_mat [measure fit $compare $reference]
# do the alignment
$all move $trans_mat

# compute the RMSD
#loop through all residues
foreach r $res {
   set ref [atomselect $mol "protein and resid $r and noh" frame 0]
   set comp [atomselect $mol "protein and resid $r and noh" frame $frame]
   set rmsd($r) [expr $rmsd($r) + [measure rmsd $comp $ref]]
   $comp delete
   $ref delete
}
}

set ave 0
foreach r $res {
   set rmsd($r) [expr $rmsd($r)/$num_steps]
   # print the RMSD
   puts "RMSD of residue $r is $rmsd($r)"
   puts $fil " $r \t $rmsd($r)"
   set res_b [atomselect $mol "resid $r"]
         $res_b set user $rmsd($r)
         $res_b delete
   set ave [expr $ave + $rmsd($r)]
}

set ave [expr $ave/[llength $res]]
puts " Average rmsd per residue: $ave"
close $fil
}

大家看看程序的问题吗?还有，之前的任何过程都没出现任何的error。先谢谢各位的指点了。

回复此楼

1楼 2010-12-15 15:38:32

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

nufang19a

金虫 (正式写手)

应助: 0 (幼儿园)
金币: 833.9
散金: 103
红花: 1
帖子: 449
在线: 101.7小时
虫号: 968803
注册: 2010-03-12
性别: GG
专业: 分子模拟

是不是我跑平衡动力学时，有些参数设计的不对啊，这是我平衡动力学的参数：
## JOB DESCRIPTION                                        ##
#############################################################
# Minimization and Equilibration of
# Ubiquitin in a Water Box

#############################################################
## ADJUSTABLE PARAMETERS                                  ##
#############################################################
structure       ../common/1kim_wb.psf
coordinates       ../common/1kim_wb.pdb
set temperature 310
set outputname    1kim_wb_eq
firsttimestep    0

#############################################################
## SIMULATION PARAMETERS                                  ##
#############################################################
# Input
paraTypeCharmm    on
parameters       ../common/par_all27_prot_lipid.inp
temperature       $temperature

# Force-Field Parameters
exclude          scaled1-4
1-4scaling       1.0
cutoff             20.
switching          on
switchdist       18.
pairlistdist       22.

# Integrator Parameters
timestep          2.0  ;# 2fs/step
rigidBonds       all  ;# needed for 2fs steps
nonbondedFreq    1
fullElectFrequency  2
stepspercycle    10

# Constant Temperature Control
langevin          on ;# do langevin dynamics
langevinDamping    5    ;# damping coefficient (gamma) of 5/ps
langevinTemp       $temperature
langevinHydrogen off ;# don't couple langevin bath to hydrogens

# Periodic Boundary Conditions
cellBasisVector1 79. 0. 0.
cellBasisVector2    0. 72. 0.
cellBasisVector3    0. 0. 72.
cellOrigin          0. 0. 0.
wrapAll          on

# PME (for full-system periodic electrostatics)
PME                yes
PMEGridSizeX       82
PMEGridSizeY       75
PMEGridSizeZ       75

# Constant Pressure Control (variable volume)
useGroupPressure    yes ;# needed for rigidBonds
useFlexibleCell    no
useConstantArea    no
langevinPiston       on
langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod  100.
langevinPistonDecay 50.
langevinPistonTemp $temperature

# Output
outputName       $outputname
restartfreq       500    ;# 500steps = every 1ps
dcdfreq          250
xstFreq          250
outputEnergies    100
outputPressure    100

#############################################################
## EXTRA PARAMETERS                                     ##
#############################################################

#############################################################
## EXECUTION SCRIPT                                     ##
#############################################################
# Minimization
minimize          5000
reinitvels       $temperature
run 20000 ;# 5ps

赞一下

回复此楼

2楼2010-12-15 15:38:50

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

nufang19a

金虫 (正式写手)

应助: 0 (幼儿园)
金币: 833.9
散金: 103
红花: 1
帖子: 449
在线: 101.7小时
虫号: 968803
注册: 2010-03-12
性别: GG
专业: 分子模拟

这个是主链的随时RMSD，使用的脚本是：
set outfile [open rmsd.dat w];
set nf [molinfo top get numframes]
set frame0 [atomselect top "protein and backbone and noh" frame 0]
set sel [atomselect top "protein and backbone and noh"]
# rmsd calculation loop
for {set i 1 } {$i < $nf } { incr i } {
$sel frame $i
$sel move [measure fit $sel $frame0]
puts $outfile "$i [measure rmsd $sel $frame0]"
}
close $outfile

赞一下

回复此楼

3楼2010-12-15 18:04:55

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zyj8119

木虫 (著名写手)

模拟EPI: 10
应助: 65 (初中生)
贵宾: 0.003
金币: 915.1
散金: 1440
红花: 35
帖子: 2936
在线: 1329.4小时
虫号: 664177
注册: 2008-11-29
性别: GG
专业: 理论和计算化学

★
小木虫(金币+0.5):给个红包，谢谢回帖交流

引用回帖:

Originally posted by nufang19a at 2010-12-15 18:04:55:
这个是主链的随时RMSD，使用的脚本是：
set outfile [open rmsd.dat w];
set nf [molinfo top get numframes]
set frame0 [atomselect top "protein and b ...

你用的是什么软件？这个GUI是什么？

赞一下(1人)

回复此楼

好好学习，天天向上。