src/gromacs/selection/tests/gensphere.py

   1 #!/usr/bin/python
   2 #
   3 # This file is part of the GROMACS molecular simulation package.
   4 #
   5 # Copyright (c) 2012, by the GROMACS development team, led by
   6 # David van der Spoel, Berk Hess, Erik Lindahl, and including many
   7 # others, as listed in the AUTHORS file in the top-level source
   8 # directory and at http://www.gromacs.org.
   9 #
  10 # GROMACS is free software; you can redistribute it and/or
  11 # modify it under the terms of the GNU Lesser General Public License
  12 # as published by the Free Software Foundation; either version 2.1
  13 # of the License, or (at your option) any later version.
  14 #
  15 # GROMACS is distributed in the hope that it will be useful,
  16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18 # Lesser General Public License for more details.
  19 #
  20 # You should have received a copy of the GNU Lesser General Public
  21 # License along with GROMACS; if not, see
  22 # http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23 # Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24 #
  25 # If you want to redistribute modifications to GROMACS, please
  26 # consider that scientific software is very special. Version
  27 # control is crucial - bugs must be traceable. We will be happy to
  28 # consider code for inclusion in the official distribution, but
  29 # derived work must not be called official GROMACS. Details are found
  30 # in the README & COPYING files - if they are missing, get the
  31 # official version at http://www.gromacs.org.
  32 #
  33 # To help us fund GROMACS development, we humbly ask that you cite
  34 # the research papers on the package. Check out http://www.gromacs.org.
  35
  36 """Script for generating spherical test configurations."""
  37
  38 import math
  39 import random
  40 import sys
  41
  42 def dot(a, b):
  43     return sum(x*y for (x, y) in zip(a, b))
  44
  45 def norm(a):
  46     return math.sqrt(dot(a, a))
  47
  48 def angle(a, b):
  49     return math.degrees(math.acos(dot(a, b) / (norm(a) * norm(b))))
  50
  51 def minangle(a, list):
  52     minangle = 180
  53     minindex = -1
  54     for index, x in enumerate(list):
  55         xangle = angle(a, x)
  56         if xangle < minangle:
  57             minangle = xangle
  58             minindex = index
  59     return (minangle, minindex)
  60
  61 def get_single_vec():
  62     while True:
  63         x = random.randint(-1000, 1000) / 1000.0
  64         y = random.randint(-1000, 1000) / 1000.0
  65         z = random.randint(-1000, 1000) / 1000.0
  66         pos = (x, y, z)
  67         dist = norm(pos)
  68         if dist <= 1.0 and dist > 0.25:
  69             return pos
  70
  71 def write_gro_title(fp, title, atomcount):
  72     fp.write(title + '\n')
  73     fp.write('%5d\n' % (atomcount))
  74
  75 def write_gro_atom(fp, resnr, resname, atomname, index, x):
  76     fp.write('%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n' %
  77             (resnr, resname, atomname, index, x[0], x[1], x[2]))
  78
  79 def write_gro_box(fp, box):
  80     fp.write('%10.5f%10.5f%10.5f\n' % box)
  81
  82 random.seed(1097)
  83 center = (0, 0, 0)
  84 cutoff = 20
  85 possamples = 500
  86 negsamples = 500
  87 totsamples = 10000
  88
  89 sys.stderr.write("Generating reference points\n")
  90 refpoints = []
  91 refpoints.append((0, 0, -1))
  92 refpoints.append((-0.5, 0.6, 0.1))
  93 refpoints.append((-0.5, -0.5, 0.25))
  94 while len(refpoints) < 30:
  95     pos = get_single_vec()
  96     if pos[0] > 0 and pos[1] > pos[0] and pos[2] > 0:
  97         refpoints.append(pos)
  98
  99 sys.stderr.write("Generating test points\n")
 100 postestpoints = []
 101 negtestpoints = []
 102 hits = 0
 103 samplecount = 0
 104 while samplecount < totsamples or len(postestpoints) < possamples or len(negtestpoints) < negsamples:
 105     pos = get_single_vec()
 106     (pangle, index) = minangle(pos, refpoints)
 107     if pangle < cutoff:
 108         hits += 1
 109         if len(postestpoints) < possamples:
 110             postestpoints.append(pos)
 111     if pangle > cutoff:
 112         if len(negtestpoints) < negsamples:
 113             negtestpoints.append(pos)
 114     samplecount += 1
 115
 116 cfrac = float(hits) / samplecount
 117 errest = math.sqrt((cfrac - cfrac*cfrac) / samplecount)
 118 sys.stderr.write('Cutoff: %f angles\n' % (cutoff))
 119 sys.stderr.write('Estimated covered fraction: %f +- %f\n' % (cfrac, errest))
 120
 121 debugfp = open('debug.txt', 'w')
 122 fp = sys.stdout
 123 count = 1 + len(refpoints) + len(postestpoints) + len(negtestpoints)
 124 write_gro_title(fp, 'Spherical test case, cutoff %f, cfrac %f +- %f' %
 125         (cutoff, cfrac, errest) , count)
 126 n = 1
 127 write_gro_atom(fp, 1, 'C', 'C', n, center)
 128 n += 1
 129 for i in range(len(refpoints)):
 130     write_gro_atom(fp, 2, 'R', 'R', n, refpoints[i])
 131     n += 1
 132 for i in range(len(postestpoints)):
 133     write_gro_atom(fp, 3, 'TP', 'TP', n, postestpoints[i])
 134     x = postestpoints[i]
 135     xangle, index = minangle(x, refpoints)
 136     refx = refpoints[index]
 137     debugfp.write('%3d%8.3f%8.3f%8.3f  %4.1f  %2d%8.3f%8.3f%8.3f\n' %
 138             (n-1, x[0], x[1], x[2], xangle, index, refx[0], refx[1], refx[2]))
 139     n += 1
 140 for i in range(len(negtestpoints)):
 141     write_gro_atom(fp, 4, 'TN', 'TN', n, negtestpoints[i])
 142     n += 1
 143 write_gro_box(fp, (10, 10, 10))
 144 debugfp.close()