应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » Gromacs/Amber/NAMD » 【原创】在gromacs的gro中插入小分子
51/1返回列表
查看: 4368  |  回复: 12
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 1 个 模拟EPI ，点击这里进行查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

zyj8119

木虫 (著名写手)

[交流] 【原创】在gromacs的gro中插入小分子 已有8人参与

在gromacs的水盒子中插入小分子一直是一个难题，今天我介绍一种编程插入方法，大家又可以看到程序的威力。可以在表面活性剂的水盒子中插入CO2,这个是原始的没有插入CO2的水盒子的gro文件，gro文件来自于杨小震的教材。
http://pic.muchong.com/file.php?id=207
在VMD中显示的图画是这样的：

FORTRAN程序是这样的：
CODE:
  integer nsurf,natom0,nco2,tatom,nch4
        real xc,yc,zc,dl,blx,bly,blz
*     blx,bly,blz are the length of simulation box in x,y,z direction, respectively.
        parameter(nsurf=36,tatom=43,natom0=9180,nch4=0,nco2=1000,
     &          blx=10,bly=10,blz=15)
        character resid(natom0)*8,atomname(natom0)*3,co2name(3)*3
        integer atomid(natom0),atomsn,molsn
        real coor(natom0,6),co2(nco2*3,3),ch4(nch4,3)
       
        open(10,file='box.gro',status='old')
        do 20 i=1,natom0
        read(10,*)resid(i),atomname(i),atomid(i),coor(i,1),coor(i,2),
     &          coor(i,3)
20    continue
      close(10)

      do 25 i=1,nco2*3
          do 26 j=1,3
          co2(i,j)=0.0
26      continue
25    continue


********************* add CO2 *************************

        do 30 i=1,nco2
50        xc=RAN2(IDUM)*blx
        yc=RAN2(IDUM)*bly
        zc=RAN2(IDUM)*blz
          do 40 j=1,natom0
          dl=sqrt((xc-coor(j,1))**2+(yc-coor(j,2))**2+(zc-coor(j,3))**2)
          if(dl.lt.0.25) goto 50
40      continue

        if(i.ge.2)then
            do 60 k=1,i-1
            kp=(k-1)*3+2
            dl=sqrt((xc-co2(kp,1))**2+(yc-co2(kp,2))**2+(zc-co2(kp,3))**2)
            if(dl.lt.0.25) goto 50
60        continue
        endif
*************** add CO2 *************************
        co2((i-1)*3+2,1)=xc
        co2((i-1)*3+2,2)=yc
        co2((i-1)*3+2,3)=zc

        co2((i-1)*3+1,1)=xc-0.115
        co2((i-1)*3+1,2)=yc
        co2((i-1)*3+1,3)=zc                      

        co2((i-1)*3+3,1)=xc+0.115
        co2((i-1)*3+3,2)=yc
        co2((i-1)*3+3,3)=zc      
30    continue


*************** add CH4 *************************

        do 62 i=1,nch4
63        xc=RAN2(IDUM)*blx
        yc=RAN2(IDUM)*bly
        zc=RAN2(IDUM)*blz
          do 64 j=1,natom0
          dl=sqrt((xc-coor(j,1))**2+(yc-coor(j,2))**2+(zc-coor(j,3))**2)
          if(dl.lt.0.25) goto 63
64      continue

            do 65 kp=1,nco2*3
            dl=sqrt((xc-co2(kp,1))**2+(yc-co2(kp,2))**2+(zc-co2(kp,3))**2)
            if(dl.lt.0.25) goto 63
65        continue


        ch4(i,1)=xc
        ch4(i,2)=yc
        ch4(i,3)=zc      
62    continue



****************** add CH4 ************************


******************OUTPUT THE NEW GRO FILE ************************

      open(70,file='system.gro',access='append')
        write(70,'(A6)')'system'
        ktotal=natom0+3*nco2+nch4
        write(70,'(I8)')ktotal
      
        do 120 i=1,nch4
        molsn=i
        atomsn=i
        write(70,666)molsn,'MET','CME',atomsn,
     &             ch4(atomsn,1), ch4(atomsn,2),ch4(atomsn,3)
120   continue

      co2name(1)='OAA'
        co2name(2)='CAC'
        co2name(3)='OAB'
        do 100 i=1,nco2
          do 110 j=1,3
          molsn=i+nch4
          atomsn=(i-1)*3+j+nch4
          kjishu=(i-1)*3+j
          write(70,666)molsn,'DRG',co2name(j),atomsn,
     &               co2(kjishu,1),co2(kjishu,2),co2(kjishu,3)
110     continue
100   continue


        do 80 i=1,nsurf
          do 90 j=1,tatom
          molsn=nco2+i+nch4
          atomsn=nco2*3+(i-1)*tatom+j+nch4
          kjishu=(i-1)*tatom+j
          write(70,666)molsn,'MOR',atomname(kjishu),atomsn,coor(kjishu,1),
     &                   coor(kjishu,2),coor(kjishu,3)
90      continue
80    continue
******************OUTPUT THE NEW GRO FILE ************************

666   format(I5,A3,4X,A3,I5,3F8.3)
      write(70,'(3F10.5)')blx,bly,blz
        write(70,*)''
      close(70)

       


        end


      FUNCTION RAN2(IDUM)
  
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
        INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11    continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END      

的操作文件，然后生成新文件，gro文件变成这样的：
http://pic.muchong.com/file.php?id=208
在VMD中的显示变成：

有一部分水分子选择了VDW后，仍然还是显示细线模型，不知道是为什么？其中nsurf是表面活性剂的个数，tatom是单个的表面活性剂分子的原子个数，nCO2是要插入CO2的个数，kjishu是表示对于某种原子个数的计数，molsn是分子个数，atomsn是总体原子个数计数。

[ Last edited by zyj8119 on 2010-12-14 at 12:33 ]
回复此楼

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

MDs-Gromacs xuexijisuan

» 本帖已获得的红花（最新10朵）

xiaowu759 lucyfly

» 猜你喜欢

» 本主题相关商家推荐: (我也要在这里推广)

» 本主题相关价值贴推荐，对您同样有帮助:
好好学习，天天向上。
1楼 2010-12-14 03:09:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

moueor

金虫 (小有名气)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
都是心血，好好支持一下！
赞一下(1人) 回复此楼
9楼2010-12-29 16:25:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 13 个回答

zyj8119

木虫 (著名写手)
引用回帖:
Originally posted by zyj8119 at 2010-12-14 03:09:42:
在gromacs的水盒子中插入小分子一直是一个难题，今天我介绍一种编程插入方法，大家又可以看到程序的威力。可以在表面活性剂的水盒子中插入CO2,这个是原始的没有插入CO2的水盒子的gro文件，gro文件来自于杨小震的教 ...

大家在不同的物质的盒子中都可以尝试这种办法。
赞一下 回复此楼
好好学习，天天向上。
2楼2010-12-14 03:11:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yalefield

金虫 (文坛精英)

老汉一枚
楼主是老杨的学生？
赞一下 回复此楼
3楼2010-12-14 08:06:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
引用回帖:
Originally posted by yalefield at 2010-12-14 08:06:57:
楼主是老杨的学生？

不是，只是拿以前的程序去做一下新的GRO文件。
赞一下 回复此楼
好好学习，天天向上。
4楼2010-12-14 12:23:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 13 个回答
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 明天应该可查了！？ +5 chengyan1220 2026-08-23 5/250 2026-08-24 23:28 by 我4大白菜
[基金申请] 如果此刻你正在为国基感到焦虑，不妨来听听这首《基金之外》 +7 scalable 2026-08-24 7/350 2026-08-24 23:01 by anata1209
[基金申请] 人气不行了 +10 fansofjerry 2026-08-21 10/500 2026-08-24 21:03 by zhanghaozhu
[基金申请] 我面上完蛋了 +9 且听虎啸 2026-08-20 10/500 2026-08-24 20:08 by 苦寒来香
[基金申请] 建议基金发布提前给出明确的时间点 +13 kulium 2026-08-21 16/800 2026-08-24 16:27 by superceng
[基金申请] 估计是周四 +4 archvillain 2026-08-18 4/200 2026-08-24 13:53 by zzuzxg
[基金申请] 放榜前的不淡定 20+4 snowwithsea 2026-08-19 17/850 2026-08-24 10:20 by echo8914667
[基金申请] 科研孤儿太难了 +18 我4大白菜 2026-08-20 19/950 2026-08-24 09:47 by 凯恩广盛催化分�
[基金申请] 什么时候开奖？ +10 CrisMessi 2026-08-18 11/550 2026-08-24 06:50 by 开心的小狮子
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +4 process2012 2026-08-23 5/250 2026-08-23 19:58 by jurkat.1640
[基金申请] 今天放榜吗？ +15 布布和一二 2026-08-19 16/800 2026-08-23 09:55 by 张春生
[基金申请] 只有每年这种时候来逛逛小木虫 +24 yaoyewhu2008 2026-08-20 26/1300 2026-08-22 17:43 by kammury
[基金申请] 时间戳今天，20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 看来今天不会放榜了？ +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 应该是下周三26日公布了吧？ +4 哈哈蛤？ 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
[基金申请] 重要消息，中午系统在维护 +11 yuleib84 2026-08-18 12/600 2026-08-20 11:09 by xskun
[基金申请] 明天放榜？ +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
信息提示
关闭
请填处理意见
关闭