24小时热门版块排行榜    

查看: 2021  |  回复: 6

nufang19a

金虫 (正式写手)


[交流] 【讨论】大家看看我的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。先谢谢各位的指点了。
回复此楼

» 收录本帖的淘帖专辑推荐

分子动力学 MDs-Gromacs gromacs

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

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

nufang19a

金虫 (正式写手)


是不是我跑平衡动力学时,有些参数设计的不对啊,这是我平衡动力学的参数:
## 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

金虫 (正式写手)


这个是主链的随时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

木虫 (著名写手)


小木虫(金币+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是什么?
好好学习,天天向上。
4楼2010-12-15 18:26:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nufang19a

金虫 (正式写手)


5楼2010-12-15 18:41:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zh1987hs

金虫 (著名写手)

分子模拟新手

★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
ghcacj(金币+2):谢谢 2010-12-17 12:43:25
引用回帖:
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 ...

namd为啥不配合VMD使用呢?很方便啊
6楼2010-12-15 22:33:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)


小木虫(金币+0.5):给个红包,谢谢回帖交流
引用回帖:
Originally posted by zh1987hs at 2010-12-15 22:33:10:

namd为啥不配合VMD使用呢?很方便啊

正解,我也正在学NAMD。
好好学习,天天向上。
7楼2010-12-15 23:12:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 nufang19a 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 投票:  有多少人是今天查系统知道结果的? +16 爱看书的可乐 2026-08-26 18/900 2026-08-30 07:50 by winsaint
[论文投稿] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +3 G6APbkg8SA6w 2026-08-29 4/200 2026-08-30 07:39 by ZPa0EcMwuECS
[考博] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 gy1nBQXYQJqL 2026-08-29 5/250 2026-08-30 07:05 by ZPa0EcMwuECS
[硕博家园] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 gy1nBQXYQJqL 2026-08-29 5/250 2026-08-30 06:55 by ZPa0EcMwuECS
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +5 ASdOkHsho7FD 2026-08-28 7/350 2026-08-30 06:10 by ZPa0EcMwuECS
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +6 ASdOkHsho7FD 2026-08-28 11/550 2026-08-30 06:00 by ZPa0EcMwuECS
[教师之家] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +5 ASdOkHsho7FD 2026-08-28 8/400 2026-08-30 05:48 by ZPa0EcMwuECS
[考研] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +4 G6APbkg8SA6w 2026-08-29 5/250 2026-08-30 02:40 by ZPa0EcMwuECS
[论文投稿] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +4 gy1nBQXYQJqL 2026-08-29 4/200 2026-08-30 02:09 by ZPa0EcMwuECS
[博后之家] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +5 gy1nBQXYQJqL 2026-08-29 6/300 2026-08-30 01:47 by ZPa0EcMwuECS
[考博] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 gy1nBQXYQJqL 2026-08-29 4/200 2026-08-30 01:34 by ZPa0EcMwuECS
[基金申请] 有没有仍没收到信息的 +5 德尚中行 2026-08-27 6/300 2026-08-29 18:10 by manplx
[基金申请] 为什么到现在没收到通知? +4 tannykie 2026-08-29 4/200 2026-08-29 17:14 by lmz0216
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +9 苏东坡二世 2026-08-23 9/450 2026-08-29 14:41 by hustersqt
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +12 Kittylucky 2026-08-27 13/650 2026-08-29 00:04 by superceng
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 我不理解! +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
信息提示
请填处理意见