Distribute Points Evenly on Surface of a Sphere

From: Rajan Vatassery (rajan_at_umn.edu)
Date: Mon Oct 18 2010 - 13:57:53 CDT

Dear List,
        I am posting the contents of a python script I wrote to do the task in
the subject line. I do not know python very well, so I understand that
it may not be the most beautiful code you've seen, but it does the job.
        It is a linux-based program, and probably will not work with windows.
Feel free to modify it at your convenience -- it is licensed under the
GPL.
        It is relevant to this list because the critical step in the program
(distributing the points evenly) needs to be done using namd.

I am posting the full script here because I'm not sure how long my
server will be active.

FULL CODE:
####################################################################

import sys
import math
import numpy
import random
import itertools

RandPts = sys.argv[1] # File to place random points into
ConfFile = sys.argv[2] # Config file to write
PSFfile = sys.argv[3] # PSF to write
numpts = int(sys.argv[4])
#numpts = int(raw_input('How Many Points? '))

f = open(RandPts,"w")
g = open(ConfFile,"w")

li = 0

outputName = 'SphDistOut'+str(numpts)+'.pdb'
paramFile = 'SphDistrNAMD.prm'
FixPts = 'SphDistrFix.pdb'

x = []
y = []
z = []

r = 2.5
Welldepth = 9999.99

#################################################################
### PDB file written, containing fixed atom at origin ###
#################################################################

# 1.00 in beta column of line below informs namd to regard CAA as fixed

f.write('REMARK Generated by $python SphDistrNAMD.py fil1 fil2 fil3 X
\n')
f.write('HETATM 1 CAA TTHP 0 0.000 0.000 0.000 0.00
1.00 C\n')

# the line below fills names array with atom names from CAA to CZZ
names = ["".join(xABC) for xABC in
itertools.product("ABCDEFGHIJKLMNOPQRSTUVWXYZ", repeat=2)]

for li in range(0,numpts):
    
    # assign random cartesian points to arrays x, y, and z
    xhld = random.random() * 2 * math.pi
    x.append(r * math.sin(xhld) * math.cos(xhld))
    yhld = random.random() * 2 * math.pi
    y.append(r * math.sin(yhld) * math.cos(yhld))
    zhld = random.random() * 2 * math.pi
    z.append(r * math.cos(zhld))

    f.write('HETATM %3d C%s TTHP 0 %- 2.3f %- 2.3f %- 2.3f
0.00 0.00 C\n' % (li+2, names[li+1], x[li], y[li], z[li]))

    li = li + 1

f.write('END\n')

f.close()

#################################################################
# Write conf file for minimization, bare minimum params given ###
#################################################################

g.write('#############################################################
\n')
g.write('#### Conf File Generated by SphereDist.py #####
\n')
g.write('#############################################################
\n')
g.write('#### Input Files Below #####
\n \n')
g.write('structure %s \n' % PSFfile)
g.write('coordinates %s \n' % RandPts)
g.write('outputName %s \n \n' % outputName)
g.write('temperature 0 \n')
g.write('#############################################################
\n')
g.write('#### Simulation Parameters #####
\n')
g.write('#############################################################
\n \n')
g.write('paraTypeCharmm on \n')
g.write('parameters %s \n\n\n' % paramFile)
g.write('fixedAtoms on \n')
# fixedAtomsCol B tells namd to look in the beta column of pdb
for
# fixed atoms mentioned in fixedAtoms on line
g.write('fixedAtomsCol B \n\n')
g.write('# Force-Field Parameters \n\n')
# must exclude none from calc, or electrostatic E will = 0
g.write('exclude none \n')
g.write('1-4scaling 1.0 \n')
g.write('cutoff 99 \n')
g.write('switching off \n')
g.write('pairlistdist 100 \n')
# binaryoutput = no will print minimized points to <outputname>.pdb.coor
g.write('binaryoutput no \n\n')
g.write('#############################################################
\n')
g.write('#### Execution Script #####
\n')
g.write('#############################################################
\n \n')
# change number below to make simulation shorter or longer as needed
g.write('minimize 10000 \n')
g.write('run 0 \n')

g.close()

#################################################################
#### Write psf file #####
#################################################################

numptplone = numpts + 1

i = open(PSFfile,"w")
i.write('PSF\n\n 3 !NTITLE\n')
i.write(' REMARKS PSF generated by SphDistrNAMD.py\n')
i.write(' REMARKS topology no_topology_file\n')
i.write(' REMARKS segment TTH { first NONE; last NONE; auto none }\n')
i.write('\n')
i.write(' %3d !NATOM \n' % numptplone)

typ = 'CG999'
# mass/charge are not related to anything physical

########################################################################
## WARNING: you need to check each of your structures so that the
## atoms are always on the surface of the sphere. If the charges are
## too great, the electrostatic component of the potential energy
surface
## will overwhelm the bond energy's component and NAMD will incorrectly
## make a second sphere outside of the first sphere during the course of
the
## minimization. This will be obvious if you look at your minimized
## structure and you notice its shape has two concentric
spheres. :WARNING
########################################################################

charge = 0.100000
mass = 15.9994
li = 1

i.write(' 1 TTH 0 TTHP CAA CG000 3.000000 15.9994
C\n')

for li in range(1,numpts+1):
    i.write(' %3d TTH 0 TTHP C%s %s % 1.6f % 2.4f
C\n' % (li + 1, names[li+1], typ, charge, mass))
    li = li + 1

# Bonds were formerly problematic. If you are having problems with your
# simulation, this is the first place to look.

flor = int(math.floor(float(numpts)/float(4)))

i.write('\n')
i.write(' %3d !NBOND: bonds\n' % numpts)
for li in range(0,flor):
    i.write(' 1 %3d 1 %3d 1 %3d
1 %3d\n' % (2 + (li*4), 3 + (li*4), 4 + (li*4), 5 + (li*4)))

if numpts%4 == 3:
    i.write(' 1 %3d 1 %3d 1 %3d\n' %
(2 + (flor*4), 3 + (flor*4), 4 + (flor*4)))
elif numpts%4 == 2:
    i.write(' 1 %3d 1 %3d\n' % (2 + (flor*4), 3
+ (flor*4)))
elif numpts%4 == 1:
    i.write(' 1 %3d\n' % (2 + (flor*4)))

i.write('\n')
i.write(' 0 !NTHETA: angles \n')
i.write('\n\n')
i.write(' 0 !NPHI: dihedrals \n')
i.write('\n\n')
i.write(' 0 !NIMPHI: impropers \n')
i.write('\n\n\n')
i.write(' 0 !NDON: donors \n')
i.write('\n\n')
i.write(' 0 !NACC: acceptors \n')
i.write('\n\n')
i.write(' 0 !NNB \n \n')

i.close()

#################################################################
##### Write Param File #####
#################################################################

ParF = open(paramFile, "w")

ParF.write('*Nothing even close to CHarmM here, just a bunch of fake
params << \n')
ParF.write('*>>>>> written by SphDistrNAMD.py for the purpose of
making<<<<<<< \n')
ParF.write('*>>>>>>>>>>>>>>>> simulation accurate as possible
<<<<<<<<<<<<<<<< \n')
ParF.write('\n\n\n')
ParF.write('BONDS \n')
ParF.write('CG000 CG999 % 7.2f % 1.4f ! from SphDistrNAMD.py \n' %
(Welldepth,r))
ParF.write('\n \n \n')
ParF.write('NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch
- \n')
ParF.write('cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin
1.5 \n')
ParF.write(' !adm jr., 5/08/91, suggested cutoff scheme
\n')
ParF.write('!see mass list above for better description of atom types
\n')
ParF.write('!hydrogens \n')
# These parameters will ensure that namd ignores any VDW contribution
ParF.write('CG999 0.0 0.0000 1.3400 ! \n')
ParF.write('CG000 0.0 0.0000 1.3400 ! \n')
ParF.write('\n \n \n')
ParF.write('END \n')

ParF.close()

##################################################################
END OF CODE

USAGE:
you will need a python interpreter to create the files to be
subsequently run in NAMD. the python command i used is:

python SphDistrNAMD.py SphereFile.pdb SphDistNAMD.conf SphDistrNAMD.psf
X (all in one line)

This will run the python routine named SphDistrNAMD.py, creating a set
of random points in SphereFile.pdb, a config file with bare minimum
parameters "SphDistNAMD.conf," and a PSF "SphDistrNAMD.psf" that will
contain (along with the pdb file) X points that represent charged
particles.

When you run the namd command on the config file, you will see that it
minimizes the coulombic repulsion between each point on the surface of a
sphere. The spherical shape is enforced by a bond between a fixed,
central atom at the origin and all the other points in the minimization.
Everything that logically can be changed in such a simulation can be
changed in the python script pasted above. I have added comments in the
script where I felt they may be helpful, so please take a look if you
have questions.

If you are not familiar with NAMD, the output of the calculation will be
in the .coor file produced after the minimization has completed.

Hope this is helpful.

Rajan

PS : also see this site for more information on this subject:
http://www.tc.umn.edu/~rajan/PtsOnSphere.html

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:56:14 CST