[交流] 【求助】求GROMACS模拟表面活性剂在油水表面的行为的例子已有3人参与

Originally posted by zyj8119 at 2010-04-02 22:25:48:
我师兄就是做这个的，不过在GROAMCS跑好的文件中，还通过自己编程，读入了吸附质的坐标。。

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=3456,tatom=1,natom0=3456,nch4=0,nco2=100,
     &          blx=4.04680,bly=4.04590,blz=4.02660)
        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='sorbent.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