Matematické Fórum

Dobrý den, snažím se naprogramovat skript v pythonu, který by měl jako výstup dát souřadnice molekul na povrchu ze souboru souřadnic molekul tvořících krychlový box vody nad kterým je vauum. Je to modifikace programu, který již fungoval pro kapku vody. Tam však nebylo nutné zavádět periodické okrajové podmínky oproti případu s boxem vody. Mám problém, že modifikovaný program stále nefunguje správně, dává moc vysoké hodnoty hustoty ( grid [i][j][k] ), nevím, kde může bát problém - zda s algoritmem zavádějícím per. ok. podm. ...  Počítají se gaussiány kolem jednotlivých grid pointů až do určitého cutoffu.
Děkuji za pomoc.
Program je zde:

import sys
import math
import numpy as np
import commands
import time

"""
Program to assign Gaussians to oxygen positions for a smooth surface definition
/DS:2013-03-25
"""

def main():   
   
    skipsymbols = np.array(['#', '@'])
    dtype       = np.float64
    name        = sys.argv[1]
   
    cmd         = 'ls ' + name
    fail,out    = commands.getstatusoutput(cmd)
    assert fail == 0
   
    f           = open(name, 'r')
   
    # read data from gromacs .gro file
    atomid, data, numatoms, xbox, ybox, zbox = read_gro(name)
    xcoord, ycoord, zcoord                   = data[0],data[1],data[2]    
   
    # define 3D grid
    Ngrid = 20  #100   # 50
   
    dx    = xbox / Ngrid
    dy    = ybox / Ngrid
    dz    = zbox / (2.0 * Ngrid)     # make this more general!
   
    grid  = [[[0.0 for k in xrange(Ngrid)] for j in xrange(Ngrid)] for i in xrange(Ngrid)] 
        grid_pt_name = [' N   '] * int(Ngrid * Ngrid * Ngrid)
   
    x    = [0.0] * int(Ngrid * Ngrid * Ngrid)
    y    = [0.0] * int(Ngrid * Ngrid * Ngrid)
    z    = [0.0] * int(Ngrid * Ngrid * Ngrid)
   
    # define parameters for density field
    sigma      = 0.21               # width of Gaussians around atom position [nm]
    nsig       = 3                  # range to add Gaussian contributions at grid point
    width      = nsig * sigma   
        width2     = width * width
    rhol       = 1.0                # liquid density in bulk [g/cm^3]
    threshdens = rhol / 2.0         # threshold for surface definition   
   
    # ... and frequently used constants
    normconst = 1.0 / ((sigma * np.sqrt(2.0 * np.pi))**3)
    tssq      = 2.0 * sigma**2
   
    # loop over oxygen coodinates
    for l in range(0, len(zcoord) - 1, 4):
       
        # loop over relevant grid points and calculate density at each grid point (add Gaussian contributions)
        for i in range(0, Ngrid - 1, 1):           
           
            distx  = ( xcoord[l] - dx * i)**2
            dist2x = ( xcoord[l] - dx * i + xbox)**2
            dist3x = (-xcoord[l] + dx * i + xbox)**2                   
            #print distx," ", dist2x," ", dist3x
            #sys.exit()       
            if ((distx > width2) and (dist2x > width2) and (dist3x > width2)):
                expx = 1000.0
            else:
                expx = 0.0
                if (distx < width2):
                    expx += distx
                if (dist2x < width2):
                    expx += dist2x
                if (dist3x < width2):   
                    expx += dist3x           
           
            #print "expx = ", expx           
           
            for j in range(0, Ngrid - 1, 1):   
               
                disty  = ( ycoord[l] - dy * j)**2
                dist2y = ( ycoord[l] - dy * j + ybox)**2
                dist3y = (-ycoord[l] + dy * j + ybox)**2   
               
                if ((disty > width2) and (dist2y > width2) and (dist3y > width2)):
                    expxy= 1000.0
                else:
                    expxy = expx
                    if (disty < width2):
                        expxy += disty
                    if (dist2y < width2):
                        expxy += dist2y
                    if (dist3y < width2):   
                        expxy += dist3y
                       
                for k in range(0, Ngrid - 1, 1):           
                   
                    distz  = ( zcoord[l] - dz * k)**2
                    dist2z = ( zcoord[l] - dz * k + zbox)**2
                    dist3z = (-zcoord[l] + dz * k + zbox)**2   
                    #print distz, "  ",dist2z, "  ", dist3z, "  ", width2
                   
                    if ((distz > width2) and (dist2z > width2) and (dist3z > width2)):
                        expxyz = 1000.0
                    else:
                        expxyz = expxy
                        if (distz < width2):
                            expxyz += distz
                        if (dist2z < width2):
                            expxyz += dist2z
                        if (dist3z < width2):   
                            expxyz += dist3z
                       
                    """                   
                    if (distz < width2):
                        grid[i][j][k] += normconst * np.exp( -(expxy + distz) / tssq)
                    if (dist2z < width2):
                        grid[i][j][k] += normconst * np.exp( -(expxy + dist2z) / tssq)
                    if (dist3z < width2):   
                        grid[i][j][k] += normconst * np.exp( -(expxy + dist3z) / tssq)
                    """
                    grid[i][j][k] += normconst * np.exp( - expxyz / tssq)
                   
                    #print " expxyz = ", expxyz
                    #print grid[i][j][k]  #, " ", i , j, k
    """   
    sys.exit()
    for i in range(0, Ngrid , 1):
        for j in range(0, Ngrid, 1):
            for k in range(0, Ngrid, 1):
                if (grid[i][j][k] > threshdens):                   
                print "N  ", 10*dx*i,"  ", 10*dy*j,"  ", 10*dz*k, "  ",grid[i][j][k]
    """
    #sys.exit()
   
    ndx_cnt = 0  #counter
   
    for i in range(0, Ngrid - 1, 1):
       
        for j in range(0, Ngrid - 1, 1):
           
            for k in range(0, Ngrid - 1, 1):
               
                if(((grid[i][j][k-1] <= threshdens) and (grid[i][j][k] >= threshdens)) or((grid[i][j][k-1] > threshdens) and (grid[i][j][k] < threshdens))):
                    zsurf = (k-1)*dz + (threshdens - grid[i][j][k-1])*dz/(grid[i][j][k]-grid[i][j][k-1])
                    x[ndx_cnt], y[ndx_cnt], z[ndx_cnt]  = i*dx, j*dy, zsurf
                    ndx_cnt += 1
                    #print "N   ", i*dx, j*dy, zsurf   # output in cartesian coods with origin one corner
                    #rad = np.sqrt((i*dx-x0)**2 + (j*dy-y0)**2 + (zsurf-z0)**2)
                    #print rad, np.arccos((zsurf-z0)/rad), np.arctan2((j*dy-y0), (i*dx-x0))
                    #print np.sqrt((i*dx-xcoord[0])**2 + (j*dy-ycoord[0])**2 + (k*dz-zcoord[0])**2)
                    #print  i*dx, "  ", j*dy, "  ", zsurf, "  ",  grid[i][j][k]
                   
                if(((grid[i][j-1][k] <= threshdens) and (grid[i][j][k] >= threshdens)) or ((grid[i][j-1][k] > threshdens) and (grid[i][j][k] < threshdens))):
                    ysurf = (j-1)*dy + (threshdens - grid[i][j-1][k])*dy/(grid[i][j][k]-grid[i][j-1][k])
                    x[ndx_cnt], y[ndx_cnt], z[ndx_cnt]  = i*dx, ysurf, k*dz
                    ndx_cnt += 1
                    #print "N   ", i * dx, ysurf, k*dz   # output in cartesian coods
                    #rad = np.sqrt((i*dx-x0)**2 + (ysurf-y0)**2 + (k*dz-z0)**2)
                    #print rad, np.arccos((k*dz-z0)/rad), np.arctan2((ysurf-y0),(i*dx-x0))
                    #print np.sqrt((i*dx-xcoord[0])**2 + (j*dy-ycoord[0])**2 + (k*dz-zcoord[0])**2)
                    #print i*dx, "  ", ysurf, "  ", k*dz, "  ",  grid[i][j][k]
                   
                if(((grid[i-1][j][k] <= threshdens) and (grid[i][j][k] >= threshdens)) or ((grid[i-1][j][k] > threshdens) and (grid[i][j][k] < threshdens))) :
                    xsurf = (i - 1) * dx + (threshdens - grid[i - 1][j][k]) * dx / (grid[i][j][k] - grid[i - 1][j][k])
                    x[ndx_cnt], y[ndx_cnt], z[ndx_cnt]  = xsurf, j*dy, k*dz
                    ndx_cnt += 1
                    #print "N   ", xsurf, j*dy, k*dz  # output in cartesian coods
                    #rad = np.sqrt((xsurf-x0)**2 + (j*dy-y0)**2 + (k*dz-z0)**2)
                    #print rad, np.arccos((k*dz-z0)/rad), np.arctan2((j*dy-y0),(xsurf-x0))
                    #print np.sqrt((i*dx-xcoord[0])**2 + (j*dy-ycoord[0])**2 + (k*dz-zcoord[0])**2)
                    #print  xsurf, "  ", j*dy, "  ", k*dz, "  ",  grid[i][j][k]
                   
    for i in range(0, Ngrid, 1):
        for j in range(0, Ngrid, 1):
            for k in range(0, Ngrid, 1):
                if (grid[i][j][k] < 0.5):
                    print "N  ", 10*i*dx, " ",10*j*dy," ", 10*k*dz
    sys.exit()
    write_gro(Ngrid*Ngrid*Ngrid, xbox, ybox, zbox, grid_pt_name, x, y, z)
   
    sys.exit()
   

####################  FUNCTION DEFINITIONS  ###########################

# reads gromacs configuration and writes lines into array
def read_gro(filename, separator=' ', skipsymbol='#'):
   
    natoms = 0   # needs to be initialized before reading in l==1
    cast = np.cast
    data = [[] for dummy in range(9)]
    l=0
   
    # reading the data from gromacs configuration file
    #print "  Reading .gro-file..."
    for line in open(filename, 'r'):
       
        bRead = True
       
        if (l==0):
            bRead = False
           
        if (l==1):
            fields = line.strip().split()
            natoms =  cast[np.int64](fields[0])
            bRead = False
           
        if (l==(natoms + 2)):
            fields = line.strip().split()
            xbox = cast[np.float64](fields[0])
            ybox = cast[np.float64](fields[1])
            zbox = cast[np.float64](fields[2])
            bRead = False
           
        if bRead:
            fields = line.strip().split()
            for i, number in enumerate(fields):
                data[i].append(number)
               
        l+=1
       
    # take care of variable types
    for i in range(0,1):   
        data[i] = cast[str](data[i])
       
    for i in range(2,8):
        data[i] = cast[np.float64](data[i])

    atomid = [0]*len(data[2])
    dat = [[0]*len(data[2]) for dummy in range(3)]
   
        ## taking care of merging columns in the .gro file
    j = 3
    for l in range(0,len(data[2])):
        if (l==9999):
            j = 2
            #print "  read_gro(): WARNING! Varying number of columns in input, control your output! Continue?"
            #time.sleep(5)
        if (l==99999):
            j = 3
            #print "  read_gro(): WARNING! Varying number of columns in input, control your output! Continue?"
            #time.sleep(5)
        if (l==109999):
            j = 2
            #print "  read_gro(): WARNING! Varying number of columns in input, control your output! Continue?"
            #time.sleep(5)
        for k in range(0,3):
            dat[k][l] = cast[np.float64](data[k+j][l])
        
        atomid[l] = cast[str](data[1][l])
                #id[l] = data[1][l]
           
    return atomid, dat, natoms, xbox, ybox, zbox
   
# reads xyz data and writes them to arrays
def read_xyz(filename, separator=' ', skipsymbol='#'):
    """ Read a xyz-file of format
        [numer of atoms]
        [Boxlenght or other comment (ignored)]
        [id  x  y  z]
    """
    cast = np.cast
    data = [[] for dummy in range(4)]
    l = 0
    for line in open(filename, 'r'):
        bRead = True
        if (l==0):
            fields = line.strip().split()
            natoms = fields[0]
            bRead = False
        if (l==1):
            fields = line.strip().split()
            xboxlen = fields[1]
            yboxlen = fields[2]
            zboxlen = fields[3]
            bRead = False
        if bRead:
            fields = line.strip().split()
            for i, number in enumerate(fields):
                data[i].append(number)
        l+=1
       
    data[0] = cast[str](data[0])
    for i in range(1,3):
        data[i] = cast[np.float64](data[i])

    return np.array(data),xboxlen,yboxlen,zboxlen,int(natoms)

# for xyz data with OHHOHH... structure, the virtual center for tip4p/2005
# is calculated and added to the arrays: OHHXOHHX...
def add_virtsites(atom_id, xpos, ypos, zpos,natoms,xbox,ybox,zbox):
    """ This routine adds the virtual sites for TIP4P/2005
        to an xyz-structure
        (Model described in gromacs topol.top files)
    """
    c = 0.13193828
    dum_id = ['MW4']*int(natoms / 3)
    xdum = [None] * int(natoms / 3)
    ydum = [None] * int(natoms / 3)
    zdum = [None] * int(natoms / 3)
    id   = [None] * int(natoms + natoms / 3)
    x    = [None] * int(natoms + natoms / 3)
    y    = [None] * int(natoms + natoms / 3)
    z    = [None] * int(natoms + natoms / 3)

    for i in range(int(natoms/3)):
        xdum[i] = float(xpos[3 * i]) + c * ( float(xpos[3 * i + 1]) + float(xpos[3 * i + 2]) - 2.0 * float(xpos[3 * i]) )
        ydum[i] = float(ypos[3 * i]) + c * ( float(ypos[3 * i + 1]) + float(ypos[3 * i + 2]) - 2.0 * float(ypos[3 * i]) )
        zdum[i] = float(zpos[3 * i]) + c * ( float(zpos[3 * i + 1]) + float(zpos[3 * i + 2]) - 2.0 * float(zpos[3 * i]) )

    # convert atom ids
    for i in range(int(natoms/3)):
        atom_id[3*i]   = 'OW1'
        atom_id[3*i+1] = 'HW2'
        atom_id[3*i+2] = 'HW3'

    # merge lists molecule-wise
    for i in range(int(natoms/3)):
        id[4*i] =  atom_id[3*i]
        id[4*i+1] = atom_id[3*i+1]
        id[4*i+2] = atom_id[3*i+2]
        id[4*i+3] = dum_id[i]
        x[4*i] =  float(xpos[3*i])/10
        x[4*i+1] = float(xpos[3*i+1])/10
        x[4*i+2] = float(xpos[3*i+2])/10
        x[4*i+3] = float(xdum[i])/10
        y[4*i] =  float(ypos[3*i])/10
        y[4*i+1] = float(ypos[3*i+1])/10
        y[4*i+2] = float(ypos[3*i+2])/10
        y[4*i+3] = float(ydum[i])/10
        z[4*i] =  float(zpos[3*i])/10
        z[4*i+1] = float(zpos[3*i+1])/10
        z[4*i+2] = float(zpos[3*i+2])/10
        z[4*i+3] = float(zdum[i])/10

    N = natoms + int(natoms/3)
   
    return id,x,y,z,N

# writes xyz coordinates to xyz file
def write_xyz(natoms,xbox,ybox,zbox,atomid,x, y, z):
    fout = open('cluster.xyz', 'w')

    # write number of atoms
    fout.write(str(natoms))
    fout.write('\n')
   
    # write box dimensions (comment in xyz)
    fout.write('BOX:'+ '  ')
    fout.write(str(xbox)+'  ')
    fout.write(str(ybox)+'  ')
    fout.write(str(zbox))
    fout.write('\n')
   
    # write atoms: id  x  y  z
    for i in range(natoms):
        fout.write(str(atomid[i])+'  ')
        fout.write(str(x[i])+'  ')
        fout.write(str(y[i])+'  ')
        fout.write(str(z[i]))
        if (i < natoms):
            fout.write('\n')
       
    fout.close()
    #return 'virtual centers attached'

# writes data to a conf.gro in proper gromacs format
def write_gro(natoms, xbox, ybox, zbox, atomid, x, y, z):
    #convert from Angstrom to nanometer
    #xbox = float(xbox)/10
    #ybox = float(ybox)/10
    #zbox = float(zbox)/10
    xbox = float(xbox)
    ybox = float(ybox)
    zbox = float(zbox)
   
    fout = open('conf.gro', 'w')
   
    # write model
    fout.write('water tip4p/2005')
    fout.write('\n')
   
    # write number of atoms
    fout.write("%5d"%(natoms))
    fout.write('\n')
   
    # write atoms: id  x  y  z
    for i in range(0,natoms):
        #fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n"%(i+1,atomid[i],i+1,x[i],y[i],z[i],0, 0, 0))
        if (((float(x[i]) <> 0) or (float(y[i]) <> 0)) or (float(z[i]) <> 0)):       
            fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f\n"%(i+1,"O",int(i),float(x[i]),float(y[i]),float(z[i])))
        #fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f\n"%(i+1,"HW2",int(4*i+1+1),float(x[4*i+1]),float(y[4*i+1]),float(z[4*i+1])))
        #fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f\n"%(i+1,"HW3",int(4*i+1+2),float(x[4*i+2]),float(y[4*i+2]),float(z[4*i+2])))
        #fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f\n"%(i+1,"MW4",int(4*i+1+3),float(x[4*i+3]),float(y[4*i+3]),float(z[4*i+3])))
        #fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f\n"%(int(i+1),str(id[i]),int(i),float(x[i]),float(y[i]),float(z[i])))
            #fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n"%(i+1,str(id[4*i+1]),int(4*i+1),float(x[4*i+1]),float(y[4*i+1]),float(z[4*i+1]),0,0,0))
            #fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n"%(i+1,str(id[4*i+2]),int(4*i+2),float(x[4*i+2]),float(y[4*i+2]),float(z[4*i+2]),0,0,0))
            #fout.write("%5dwater%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n"%(i+1,str(id[4*i+3]),int(4*i+3),float(x[4*i+3]),float(y[4*i+3]),float(z[4*i+3]),0,0,0))
       
    # write box dimensions (comment in xyz)
    fout.write("%10.5f%10.5f%10.5f\n"%(float(xbox),float(ybox),float(zbox)))
   

    fout.close()
   
    print '  wrote .gro file.'
    return '  wrote .gro file.'

main()