24小时热门版块排行榜    

查看: 2027  |  回复: 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的回帖

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的回帖
查看全部 7 个回答

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的回帖
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 国社科又开始会评了,不知道这次命运如何 +7 雨打竹帘 2026-08-30 11/550 2026-08-31 23:16 by hittle2008
[基金申请] 投票:  有多少人是今天查系统知道结果的? +17 爱看书的可乐 2026-08-26 19/950 2026-08-31 21:30 by xiangy672
[基金申请] 基金未中,这种答复是模板吗? +6 zhaosm1982 2026-08-27 7/350 2026-08-31 21:18 by qdxxmc
[基金申请] 哪位高人中了,把查询到的截图贴出来让我看看,让我长长见识 +6 yuleib84 2026-08-26 7/350 2026-08-31 19:46 by 鱼翔浅底1
[基金申请] 能否申诉? +7 echo8914667 2026-08-30 8/400 2026-08-31 17:00 by yihongxu
[文学芳草园] 梦想 +5 myrtle 2026-08-26 5/250 2026-08-31 14:40 by hahaboy
[基金申请] 中青基了要发朋友圈吗? +7 349506619 2026-08-28 7/350 2026-08-31 13:39 by 冼亮淀粉酶
[基金申请] 29号明天会评吗 +4 笨笨唐 2026-08-28 4/200 2026-08-31 09:30 by huixian257
[基金申请] 有没有仍没收到信息的 +7 德尚中行 2026-08-27 8/400 2026-08-30 20:52 by purplejack
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +15 晴天加油 2026-08-26 16/800 2026-08-29 18:28 by symmetry
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 看板上这么多中的,有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 我不理解! +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 某些机构,以效率低为荣,以效率低作为存在感 +9 yuleib84 2026-08-25 10/500 2026-08-25 17:14 by alexon
信息提示
请填处理意见