24小时热门版块排行榜    

查看: 2022  |  回复: 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

金虫 (正式写手)


这个是主链的随时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的回帖
查看全部 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的回帖

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的回帖

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的回帖
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +3 G6APbkg8SA6w 2026-08-29 4/200 2026-08-30 07:28 by ZPa0EcMwuECS
[考研] 售SCI一区T0P文章,我:8.O.55.1.O.54,科目齐全,可+急 +3 gy1nBQXYQJqL 2026-08-29 4/200 2026-08-30 06:57 by ZPa0EcMwuECS
[考博] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +6 ASdOkHsho7FD 2026-08-28 10/500 2026-08-30 06:11 by ZPa0EcMwuECS
[公派出国] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +4 gy1nBQXYQJqL 2026-08-29 4/200 2026-08-30 01:45 by ZPa0EcMwuECS
[考博] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 gy1nBQXYQJqL 2026-08-29 4/200 2026-08-30 01:34 by ZPa0EcMwuECS
[公派出国] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +5 ASdOkHsho7FD 2026-08-28 7/350 2026-08-30 01:12 by ZPa0EcMwuECS
[论文投稿] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +4 ASdOkHsho7FD 2026-08-28 5/250 2026-08-30 00:41 by ZPa0EcMwuECS
[基金申请] 为什么到现在没收到通知? +4 tannykie 2026-08-29 4/200 2026-08-29 17:14 by lmz0216
[基金申请] 系统查不到 +11 董八千 2026-08-26 11/550 2026-08-28 18:06 by Leogzhya
[基金申请] 基金不中,共勉 +11 eulota 2026-08-26 11/550 2026-08-28 14:22 by 火星超人xi
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 哪位高人中了,把查询到的截图贴出来让我看看,让我长长见识 +5 yuleib84 2026-08-26 6/300 2026-08-28 00:02 by yudaoqian88
[基金申请] 怎么看青基中了没有啊 +5 叶九微 2026-08-26 5/250 2026-08-27 10:35 by l_zh2008
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 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
[基金申请] 国际合作可查了,中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
[基金申请] 明天应该可查了!? +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 如果此刻你正在为国基感到焦虑,不妨来听听这首《基金之外》 +8 scalable 2026-08-24 8/400 2026-08-25 12:52 by jnhyjjm
信息提示
请填处理意见