24小时热门版块排行榜    

CyRhmU.jpeg
查看: 1441  |  回复: 5
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

wangcui2011

金虫 (小有名气)

[求助] NWChem计算生成的.top文件和 .rst文件各行各列的含义

各位前辈,Topology 文件和 Restart 文件是NWChem做QMMM运算需要的两个很重要的文件,想问下,这两种文件里面各行各列都是表示什么含义。如果哪位虫友知道,希望可以帮忙解答。先说一声:谢谢!
我下面给出一个例子,这个例子是nwchem documentation里面提供的:
(1) cloh.top
NWChem topology file                                                            
                                                                                
Generated by the NWChem prepare module                                          
    4.60000001/28/09   22:55:38 amber                       
    1
    5
    4
    2
    0.833333    1.000000    0.000000    0
    1   17 CL   Q   35.453000
    2    8 OH   Q   15.999400
    3    1 HO   Q    1.008000
    4    8 OWS  w   15.999400
    5    1 HWS  w    1.008000
    1    1 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    1    2 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    1    3 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    1    4 6.11777E-03 3.05888E-03 1.79398E-05 8.96989E-06
    1    5 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    2    2 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    2    3 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    2    4 2.71815E-03 1.35907E-03 2.83540E-06 1.41770E-06
    2    5 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    3    3 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    3    4 4.31389E-04 2.15694E-04 1.90177E-07 9.50885E-08
    3    5 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    4    4 2.61691E-03 1.30845E-03 2.63324E-06 1.31662E-06
    4    5 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    5    5 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    1    0.000000 0.00000E+00
    2   -0.847600 0.00000E+00
    3    0.423800 0.00000E+00
    4    0.000000 0.00000E+00
      1CLD             1   -2      1      1    1    1
      2OH              2   -2      2      2    1    1
      3      3      0      0      0         0         0
      3      0      0      0      0         0         0         0    0    0
spce       OW                                4    2    2    2      1
spce      2HW                                5    3    3    3      2
spce      3HW                                5    3    3    3      3
      1      2      1      1
    0.100000 1.00000E+06
      1      3      1      2
    0.100000 1.00000E+06
      2      3      1      3
    0.163333 1.00000E+06
spce      HOH
    5.220000
CLD        Cl     1      1      1      1      1    1    4    4    4    0    0  1      1
OH         O1     1      2      2      2      2    2    4    4    4    0    0  1      2
OH         H1     1      2      2      2      2    3    4    4    4    0    0  1      3

(2) cloh_rs.rst
Restart file


    4.20000001/29/09   14:08:50     1    F
                  01/29/09   14:08:50           cloh_opt00                                                  
    1    0    0.097665
    2.624300    0.000000    0.000000
    0.000000    2.624300    0.000000
    0.000000    0.000000    2.624300
0.00000E+00
    0.000000    0.000000    0.000000
       590         3         2         3       590         1         2    0    0
     0.08394658   0.18596279  -0.17665358   0.00000000   0.00000000   0.00000000
    -0.01256243   0.20438730  -0.15803778   0.00000000   0.00000000   0.00000000
     0.09813034   0.08731282  -0.18483929   0.00000000   0.00000000   0.00000000
00   0.00000000   0.00000000   0.00000000
     0.06624620  -0.10225782  -0.22140211   0.00000000   0.00000000   0.00000000
     0.06556329  -0.10534104  -0.12145195   0.00000000   0.00000000   0.00000000
    -0.02771490  -0.09757372  -0.25530455   0.00000000   0.00000000   0.00000000
..   ..........    .........  ...........    ..........  ...........  ...........


00   0.00000000   0.00000000   0.00000000
00  -0.15572466   0.02643901   0.02211421   0.00000000   0.00000000   0.00000000    0
00   0.07232217  -0.05744605   0.06064107   0.00000000   0.00000000   0.00000000    0
00   0.10481826   0.03465585   0.05711075   0.00000000   0.00000000   0.00000000    0
     0.00000000   0.00000000   0.00000000
     0.00000000   0.00000000   0.00000000
  1  1
restart input
      1      1
      1      1
      0      1      0      0
      0      0      0      0      0   1000    500
    0.000000    0.001000
    1.000000    1.000000
    100    0.000001
    100    0.000001
    0 0.10250E+06    0.500000 0.45300E-09    0
    0  298.150000    0.100000    0.100000  298.150000    0.000000    0.000000
      0      0  298.150000       12345
     10    100      0      0      0
      0      1
      0      0      0      0      0      0
      1
      0      0
      0      0      0      0   1000      0      0      0
    0.000000    0.000000
      0      0
      0      0
      1      0    0.000000
      1      0    0.000000
restart properties
    150      0      0
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00
  ..................  ...................  .................. ....................






restart space
     32     33      2      4      4      2      4      4     32
  0.131215000000E+01  0.262430000000E+01
  0.656075000000E+00  0.131215000000E+01  0.196822500000E+01  0.262430000000E+01
  0.656075000000E+00  0.131215000000E+01  0.196822500000E+01  0.262430000000E+01
      0     11
         0         0         0         0
         2         2         0       243
         4         4         0       543
         8         8         0       695
        10        10         0      1016
        12        12         0      1241
        16        16         0      1357
        18        18         0      1500
        20        20         0      1608
        26        26         0      1641
        28        28         0      1890
1     27
         1         1         1         0
         0         0         1       228
         2         2         1       472
         3         3         1       650
         4         4         1       931
         5         5         1      1040
         6         6         1      1169
         8         8         1      1361
         9         9         1      1542
        10        10         1      1886
        11        11         1      2018
        12        12         1      2209
        13        13         1      2315
        14        14         1      2430
        16        16         1      2566
        17        17         1      2675
        18        18         1      2813
        19        19         1      2915
        20        20         1      3023
        21        21         1      3056
        22        22         1      3089
        24        24         1      3192
        26        26         1      3378
        27        27         1      3527
        28        28         1      3740
        29        29         1      3844
        30        30         1      3957
      2     11
         2         2         2         0
         4         4         2       203
         6         6         2       523
        10        10         2       646
        12        12         2       944
        14        14         2      1128
        18        18         2      1231
        20        20         2      1342
        22        22         2      1452
        28        28         2      1485
        30        30         2      1674
       ...........................................
回复此楼
想一千次,不如去做一次。华丽的跌倒,胜过无谓的徘徊。千里之行,始于足下,即使无功而返,也不浪费大好青春。
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wangcui2011

金虫 (小有名气)

送鲜花一朵
引用回帖:
4楼: Originally posted by dulikai at 2012-09-18 19:42:29
#! /bin/python

# This script converts amber rst to nwchem rst format.

class NwFmtRst():
        def __init__(self, filename='crown.rst'):
                # file pointer
                if filename == '':
                        self.filename = & ...

谢谢大侠这么热心的回帖。虽然其中的很多代码看不太懂,但是大致的意思还是能明白些。相信你的这个回复对很多学习这个软件的同学也会有很大的帮助
能问一下,上面的脚本是你自己写的,还是安装程序的过程中生成的?恕我在程序方面的无知
对了,麻烦你下次回帖的时候,选择“应助回帖”,不然我没办法给你金币
想一千次,不如去做一次。华丽的跌倒,胜过无谓的徘徊。千里之行,始于足下,即使无功而返,也不浪费大好青春。
5楼2012-09-19 11:04:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答

wangcui2011

金虫 (小有名气)

Restart 文件太长,只给出部分的,不过所有的restart file大致都包含那几块。如果哪位好心的虫友知道,帮忙解答一下。谢谢了。
想一千次,不如去做一次。华丽的跌倒,胜过无谓的徘徊。千里之行,始于足下,即使无功而返,也不浪费大好青春。
2楼2012-09-17 16:55:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dulikai

铁虫 (初入文坛)

看看这个有用不
#! /bin/python

# This script converts amber rst to nwchem rst format.

class NwFmtRst():
        def __init__(self):
                self.overflow_flag = 0
                self.filename = "nw.rst"
                self.topname = "nw.top"
                self.title = ""
                self.hist = ""
                self.atomnum = -1
                self.MAXATOM = 100000
               
                self.rst = {}
                self.rst['coord'] = 1
                self.rst['vel'] = 1
                self.rst['box'] = 1
               
                self.atom = {}
                self.atom['coord'] = []
                self.atom['vel'] = []
                self.atom['box'] = []
               
                self.dim = {}
                # crystal solvent(c) and not crystal solvent(x)
                # element key: coord, vel, iwm
                self.wmc = {}
                self.wmc['coord'] = []
                self.wmc['vel'] = []
                self.wmc['iwm'] = []
                self.wmx = {}
                self.wmx['coord'] = []
                self.wmx['vel'] = []
                self.wmx['iwm'] = []
                # solute atom
                # dic key: coord, vel, isar, ips
                self.solute = {}
                self.solute['coord'] = []
                self.solute['vel'] = []
                self.solute['isar'] = []
                self.solute['ips'] = []
               
                # @ Initialize the top related variables.
                self.top = {}
                self.top['dim'] = {}
               
               
                print """
                -------------------------------------------------------
                @Comment: This is amber/nwchem Interface.
                Just a format Converter.
                @Address: Come form chem306.
                @Date: 2010/03/25
                -------------------------------------------------------
                """

        def __rd_rst_atom(self, fin):
                """
                        read a section in rst/crd file, i.e. coordinate/vel ...
                        This is a private type func...
                """
                # read amb atom coord and vel (optional)
                natom = self.atomnum
                atom = []
                i = 0
               
                line = fin.readline()
                nline = len(line)
                while line!="" and i < natom:
                        # now, a very strange test:
                        # see ./src/antechamber/rst.c for reference. &_&
                        if line[4] == ".":
                                x = float(line[0:12])
                                y = float(line[12:24])
                                z = float(line[24:36])
                                coord = [x, y, z]
                                atom.append(coord)
                                i = i+1
                        else:
                                print "format error"
                                break
                               
                                # then, again another three
                        if nline > 48 and line[40] == ".":
                                x = float(line[36:48])
                                y = float(line[48:60])
                                z = float(line[60:72])
                                coord = [x, y, z]
                                atom.append(coord)
                                i = i+1
                        else:
                                print "Maybe, end of file. [X]"
                                print "current line is :\n" + line
                                break;
                               
                        line = fin.readline()
                        nline = len(line)

                print "end of this section"
                print "total number of atom is %d, and current number is %d"  % (natom, i)
               
                return atom

        def __rd_rst_dim(self, fin):
                """
                read in dimensional info.
                """
                # @comment: dimensions on the restart file
                # c
                # c     1 i10 nwm   number of solvent molecules
                # c     2 i10 nwa   number of atoms per solvent molecule
                # c     3 i10 nsm   number of solute molecules
                # c     4 i10 nsa   number of solute atoms
                # c     5 i10 nwmc  number of crystal solvent molecules
                # c     6 i10 nsf   number of solute fractions
                # c     7 i10 nss   number of solute segments
                # c     8 i5  nprev number of processors used in previous job
                # c     9 i5  noe   number of noe constraints
                # c
                # c         6 & 8 may be useless, and rdrst do not read it
                # c
                line = fin.readline()
                dim = {}
                dim['nwm'] = int(line[0:10])
                dim['nwa'] = int(line[10:20])
                dim['nsm'] = int(line[20:30])
                dim['nsa'] = int(line[30:40])
                dim['nwmc'] = int(line[40:50])
                dim['nsf'] = int(line[50:60])
                dim['nss'] = int(line[60:70])
                dim['nprev'] = int(line[70:75])
                dim['noe'] = int(line[75:80])
               
                # print dim
               
                return dim
               
        def __rd_rst_nwm(self, fin, xc):
                """
                read in solvent molecules
                crystal or not
                """
                nwm = self.dim['nwm']
                nwmc = self.dim['nwmc']
                nwa = self.dim['nwa']
                w = {}
                w['coord'] = []
                w['vel'] = []
                w['iwm'] = []
                # checking ...
                if xc == "c":
                        nw = nwmc
                else:
                        nw = nwm - nwmc
               
                if nw <= 0:
                        print xc + " water is zero: Noting Done"
                       
                # reading...
                if nwm > 0:
                        for i in range(nw):
                                for j in range(nwa):
                                        line = fin.readline()
                               
                                        cx = float(line[2:15])
                                        cy = float(line[15:28])
                                        cz = float(line[28:41])
                                        vx = float(line[41:54])
                                        vy = float(line[54:67])
                                        vz = float(line[67:80])
                                        coord = [cx, cy, cz]
                                        vel = [vx, vy, vz]

                                        w['coord'].append(coord)
                                        w['vel'].append(vel)
                                        # if lforces, IGNORE
                               
                                # read iwrc
                                line = fin.readline()
                                i = int(line[0:1])
                                w['iwm'].append(i)
                               
                return w               
               
        def __rd_rst_solute(self, fin):
                """
                read in solute molecules
                """
                nsa = self.dim['nsa']
                s = {}
                s['coord'] = []
                s['vel'] = []
                s['isar'] = []
                s['ips'] = []
                # checking ...
                if nsa <= 0:
                        print "Number of solute atom is zero: Noting Done"
                # reading...
                if nsa > 0:
                        for i in range(nsa):
                                line = fin.readline()
                               
                                cx = float(line[2:15])
                                cy = float(line[15:28])
                                cz = float(line[28:41])
                                vx = float(line[41:54])
                                vy = float(line[54:67])
                                vz = float(line[67:80])
                               
                                isar = int(line[0:1])
                                ips = int(line[80:85])
                               
                                coord = [cx, cy, cz]
                                vel = [vx, vy, vz]
                                # if lforces, IGNORE
                               
                               
                                s['coord'].append(coord)
                                s['vel'].append(vel)
                                s['isar'].append(isar)
                                s['ips'].append(ips)
                               
                return s               



               
        def rd_rst(self, filename="nw.rst":
                """
                        open and read amber rst file.
                        this can be use for nwchem version >=4.6
                        thanks for pre_rdrst.F#rdrst and pre_wrtrst#wrtrst
                        and for version 3.1, maybe one can refer the sub rrst and wtrst.
                """
                self.filename = filename
               
                try:
                        fin = open(filename, "r"
                except:
                        print "Error: while opening " + filename
                        exit(1)
               
                ### read in nw rst file       
                # read nw title; jump three line
                line=fin.readline()               
                self.title += line
                line=fin.readline()               
                self.title += line
                line=fin.readline()               
                self.title += line
                # print self.title
               
                # read some tedious info.
                # rdummy means useless var. here
                # copy from pre_rdrst.F in prepar@nwchem
                line=fin.readline()               
                rdummy = float(line[0:12])
                rstdat = line[12:22]
                rdtim = line[22:32]
                nhist = int(line[32:37])
                lforces = line[41:42]
               
                # print rstdat, rdtim, nhist, lforces

                # U know, i have read nhist, which is the record for previous calc. type
                # then if nhist > 0, i should read the hist entity. Now checking and reading hist...
                if nhist > 0:
                        for i in range(nhist):
                                dummy = fin.readline()
                                self.hist += dummy
                # print self.hist
               
                # read in pbc type[period/vacuum], box type, and box dimension
                # i think argos md code is not in a good code style. Bad bad very bad...
                # Nwchem argos md module is tedious... pay more attention to understand it, if you want.
                line = fin.readline()
                npbtyp = int(line[0:5])
                nbxtyp = int(line[5:10])
                rsgm = float(line[10:22]) # This variable is not useful in current version. bad
               
                # read box dimen. remember to jump the useless zero
                line = fin.readline()
                x = float(line[0:12])
                line = fin.readline()
                y = float(line[12:24])
                line = fin.readline()
                z = float(line[24:36])
               
                # print x, y, z
               
                # read a dummy line, contain format(e12.5)
                line = fin.readline()
                dummy = float(line[0:12])
                # print dummy
                # read another dummy line, contain format(3f12.6)
                line = fin.readline()
               
                # @ read the dimensional info.
                self.dim = self.__rd_rst_dim(fin)
               
                # @ Finally, now, we can read the coordinates.
                # ... read crystal/not crystal solvent first
                self.nwc = self.__rd_rst_nwm(fin, "c"
                self.nwx = self.__rd_rst_nwm(fin, "x"
               
                # read solute atom
                self.solute = self.__rd_rst_solute(fin)
                # read solute molecule, IGNORE in this verion
               
                # read solute seq, IGNORE in this verion
               
                fin.close()
               
                return       

               

        def wrt_rst(self, filename="nwout.rst":
                """
                write out amber rst file. haha
                """       
                ### Exception checking...
                try:
                        fout = open(filename, "w"
                except:
                        print "Error: while opening " + filename
                        exit(1)
                       
                fout.close()
               
                return

       
                       
rst = NwFmtRst()
rst.rd_rst("crown.rst"
rst.wrt_rst()
3楼2012-09-18 19:34:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dulikai

铁虫 (初入文坛)

#! /bin/python

# This script converts amber rst to nwchem rst format.

class NwFmtRst():
        def __init__(self, filename='crown.rst'):
                # file pointer
                if filename == '':
                        self.filename = "crown.rst"
                else:
                        self.filename = filename
               
                self.fp = {'input': '', 'output': ''}
                self.fpin = self.fp['input']
                self.fpout = self.fp['output']
               
                self.atomnum = -1
                self.MAXATOM = 100000
               
                # data structure of rst file format
                self.rst = {}
                # title
                self.rst['title'] = []
                self.rst['hist'] = {}
               
                self.rst['coord'] = 1
                self.rst['vel'] = 1
                self.rst['box'] = 1
               
                self.atom = {}
                self.atom['coord'] = []
                self.atom['vel'] = []
                self.atom['box'] = []
               
                self.dim = {}
                # crystal solvent(c) and not crystal solvent(x)
                # element key: coord, vel, iwm
                self.wmc = {}
                self.wmc['coord'] = []
                self.wmc['vel'] = []
                self.wmc['iwm'] = []
                self.wmx = {}
                self.wmx['coord'] = []
                self.wmx['vel'] = []
                self.wmx['iwm'] = []
                # solute atom
                # dic key: coord, vel, isar, ips
                self.solute = {}
                self.solute['coord'] = []
                self.solute['vel'] = []
                self.solute['isar'] = []
                self.solute['ips'] = []
               
                # @ Initialize the top related variables.
                self.top = {}
                self.top['dim'] = {}
               
               
                print """
                -------------------------------------------------------
                @Comment: This is amber/nwchem Interface.
                Just a format Converter.
                @Address: Come form chem306.
                @Date: 2012/03/25
                -------------------------------------------------------
                """

        def __open_file(self, filename, mode='r'):
                """
                        ### Exception checking...
                """
                try:
                        fp = open(filename, mode)
                except:
                        print "Error: while opening " + filename
                        exit(1)
               
                return fp
               

        # >S rd/wrt title $
        def __rd_rst_title(self):
                """
                        first three line about title info.
                """
                fpin = self.fpin
               
                # read nw title; jump three line
                self.rst['title'].append(fpin.readline()[:-1])
                self.rst['title'].append(fpin.readline()[:-1])
                self.rst['title'].append(fpin.readline()[:-1])
               
                return
               
               
        def __rd_rst_hist(self):
                """
                        the history info. about previous calc.
                """
                fpin = self.fpin
               
                # read nw version, date, time, nhist, lforces, and hist...
                # rdummy means useless var. here
                #
                line=fpin.readline()               
                self.top['hist']['version'] = float(line[0:12])
                self.top['hist']['rstdate'] = line[12:22]
                self.top['hist']['rdtime'] = line[22:32]
                self.top['hist']['nhist'] = int(line[32:37])
                self.top['hist']['lforces'] = line[41:42]
                self.top['hist']['content'] = []
               
                # i have read nhist, which is the record for previous calc. type
                # then if nhist > 0, i should read the hist entity. Now checking and reading hist...
                if nhist > 0:
                        for i in range(nhist):
                                dummy = fin.readline()
                                self.top['hist']['content'].append(dummy)
               
                return

        def __rd_rst_pbc(self):
                """
                        read in pbc type[period/vacuum], box type, and box dimension
                """
                # % pbctype = npbtyp in pre_wrtrst.F  
                # c 1 for periodic, while 0 for vacuum.
                # c
                # % boxtype = nbxtyp
                # c 0 for cube/box/sphere, 2 for troct : [default: 0]
                # % rsgm
                # c see manual md section : program control options.
                # segmentsize
                # sets the maximum size of a segment. This value is used to
                # determine which segments at the boundary of the cutoff radius                
                # should be considered in the generation of the pairlists.
                # This value is also determined by the prepare module and
                # written to the restart file.
                # Use of this directive is not needed for simulations that use the current
                # prepare module to generate the restart file.               
                # % boxsize
                # box size info.
                # 2.000000    0.000000    0.000000
                # 0.000000    2.000000    0.000000
                # 0.000000    0.000000    2.000000
                # c
                # Nwchem argos md module is tedious...
                # pay more attention to understand it, if you want.
                # c
                line = fin.readline()
                pbctype = int(line[0:5])
                boxtype = int(line[5:10])
                rsgm = float(line[10:22])
               
                # read box dimen. remember to jump the useless zero
                x = float(fin.readline()[0:12])
                y = float(fin.readline()[12:24])
                z = float(fin.readline()[24:36])
               
                self.rst['sys']['pbctype'] = pbctype
                self.rst['sys']['boxtype'] = boxtype
                self.rst['sys']['rsgm'] = rsgm
                self.rst['sys']['box'] = [x, y, z]
               
                # print x, y, z
               
                return
               
               
               
        def __rd_rst_atom(self, fin):
                """
                        read a section in rst/crd file, i.e. coordinate/vel ...
                        This is a private type func...
                """
                # read amb atom coord and vel (optional)
                natom = self.atomnum
                atom = []
                i = 0
               
                line = fin.readline()
                nline = len(line)
                while line!="" and i < natom:
                        # now, a very strange test:
                        # see ./src/antechamber/rst.c for reference. &_&
                        if line[4] == ".":
                                x = float(line[0:12])
                                y = float(line[12:24])
                                z = float(line[24:36])
                                coord = [x, y, z]
                                atom.append(coord)
                                i = i+1
                        else:
                                print "format error"
                                break
                               
                                # then, again another three
                        if nline > 48 and line[40] == ".":
                                x = float(line[36:48])
                                y = float(line[48:60])
                                z = float(line[60:72])
                                coord = [x, y, z]
                                atom.append(coord)
                                i = i+1
                        else:
                                print "Maybe, end of file. [X]"
                                print "current line is :\n" + line
                                break;
                               
                        line = fin.readline()
                        nline = len(line)

                print "end of this section"
                print "total number of atom is %d, and current number is %d"  % (natom, i)
               
                return atom

        def __rd_rst_dim(self, fin):
                """
                read in dimensional info.
                """
                # @comment: dimensions on the restart file
                # c
                # c     1 i10 nwm   number of solvent molecules
                # c     2 i10 nwa   number of atoms per solvent molecule
                # c     3 i10 nsm   number of solute molecules
                # c     4 i10 nsa   number of solute atoms
                # c     5 i10 nwmc  number of crystal solvent molecules
                # c     6 i10 nsf   number of solute fractions
                # c     7 i10 nss   number of solute segments
                # c     8 i5  nprev number of processors used in previous job
                # c     9 i5  noe   number of noe constraints
                # c
                # c         6 & 8 may be useless, and rdrst do not read it
                # c
                line = fin.readline()
                dim = {}
                dim['nwm'] = int(line[0:10])
                dim['nwa'] = int(line[10:20])
                dim['nsm'] = int(line[20:30])
                dim['nsa'] = int(line[30:40])
                dim['nwmc'] = int(line[40:50])
                dim['nsf'] = int(line[50:60])
                dim['nss'] = int(line[60:70])
                dim['nprev'] = int(line[70:75])
                dim['noe'] = int(line[75:80])
               
                # print dim
               
                return dim
               
        def __rd_rst_nwm(self, fin, xc):
                """
                read in solvent molecules
                crystal or not
                """
                nwm = self.dim['nwm']
                nwmc = self.dim['nwmc']
                nwa = self.dim['nwa']
                w = {}
                w['coord'] = []
                w['vel'] = []
                w['iwm'] = []
                # checking ...
                if xc == "c":
                        nw = nwmc
                else:
                        nw = nwm - nwmc
               
                if nw <= 0:
                        print xc + " water is zero: Noting Done"
                       
                # reading...
                if nwm > 0:
                        for i in range(nw):
                                for j in range(nwa):
                                        line = fin.readline()
                               
                                        cx = float(line[2:15])
                                        cy = float(line[15:28])
                                        cz = float(line[28:41])
                                        vx = float(line[41:54])
                                        vy = float(line[54:67])
                                        vz = float(line[67:80])
                                        coord = [cx, cy, cz]
                                        vel = [vx, vy, vz]

                                        w['coord'].append(coord)
                                        w['vel'].append(vel)
                                        # if lforces, IGNORE
                               
                                # read iwrc
                                line = fin.readline()
                                i = int(line[0:1])
                                w['iwm'].append(i)
                               
                return w               
               
        def __rd_rst_solute(self, fin):
                """
                read in solute molecules
                """
                nsa = self.dim['nsa']
                s = {}
                s['coord'] = []
                s['vel'] = []
                s['isar'] = []
                s['ips'] = []
                # checking ...
                if nsa <= 0:
                        print "Number of solute atom is zero: Noting Done"
                # reading...
                if nsa > 0:
                        for i in range(nsa):
                                line = fin.readline()
                               
                                cx = float(line[2:15])
                                cy = float(line[15:28])
                                cz = float(line[28:41])
                                vx = float(line[41:54])
                                vy = float(line[54:67])
                                vz = float(line[67:80])
                               
                                isar = int(line[0:1])
                                ips = int(line[80:85])
                               
                                coord = [cx, cy, cz]
                                vel = [vx, vy, vz]
                                # if lforces, IGNORE
                               
                               
                                s['coord'].append(coord)
                                s['vel'].append(vel)
                                s['isar'].append(isar)
                                s['ips'].append(ips)
                               
                return s               



               
        def rd_rst(self):
                """
                        open and read amber rst file.
                        this can be use for nwchem version >=4.6
                        thanks for pre_rdrst.F#rdrst and pre_wrtrst#wrtrst
                        and for version 3.1, maybe one can refer the sub rrst and wtrst.
                """
               
                self.fpin = self.__open_file(self.filename, mode='r')
                ### read in nw rst file       
                # read nw title; jump three line
                self.__rd_rst_title()
                # read version, time .. hist. info.
                self.__rd_rst_hist()
               
               
               
               
                # read in pbc type[period/vacuum], box type, and box dimension
                # i think argos md code is not in a good code style. Bad bad very bad...
                # Nwchem argos md module is tedious... pay more attention to understand it, if you want.
                line = fin.readline()
                npbtyp = int(line[0:5])
                nbxtyp = int(line[5:10])
                rsgm = float(line[10:22]) # This variable is not useful in current version. bad
               
                # read box dimen. remember to jump the useless zero
                line = fin.readline()
                x = float(line[0:12])
                line = fin.readline()
                y = float(line[12:24])
                line = fin.readline()
                z = float(line[24:36])
               
                # print x, y, z
               
                # read a dummy line, contain format(e12.5)
                line = fin.readline()
                dummy = float(line[0:12])
                # print dummy
                # read another dummy line, contain format(3f12.6)
                line = fin.readline()
               
                # @ read the dimensional info.
                self.dim = self.__rd_rst_dim(fin)
               
                # @ Finally, now, we can read the coordinates.
                # ... read crystal/not crystal solvent first
                self.nwc = self.__rd_rst_nwm(fin, "c"
                self.nwx = self.__rd_rst_nwm(fin, "x"
               
                # read solute atom
                self.solute = self.__rd_rst_solute(fin)
                # read solute molecule, IGNORE in this verion
               
                # read solute seq, IGNORE in this verion
               
                fin.close()
               
                return       

               

        def wrt_rst(self, filename="nwout.rst":
                """
                write out amber rst file. haha
                """       
                ### Exception checking...
                try:
                        fout = open(filename, "w"
                except:
                        print "Error: while opening " + filename
                        exit(1)
                       
                fout.close()
               
                return

       
                       
rst = NwFmtRst("crown.rst"
rst.rd_rst()
rst.wrt_rst()

» 本帖已获得的红花(最新10朵)

4楼2012-09-18 19:42:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
信息提示
请填处理意见