[ViewVC] Annotation of: OpenMD/trunk/samples/Madelung/dipoles/buildDipolarArray

#! /usr/bin/env python

"""Dipolar Lattice Builder

Creates cubic lattices of dipoles to test the dipole-dipole
interaction code.  

Usage: buildDipolarArray 

Options:
  -h, --help              show this help
  -x, --array-type-X      use one of the basic "X" arrays
  -y, --array-type-Y      use one of the basic "Y" arrays
  -z, --array-type-Z      use one of the basic "Z" arrays
  -a, --array-type-A      use array type "A" (default)
  -b, --array-type-B      use array type "B"
  -l, --lattice=...       use the specified lattice ( SC, FCC, or BCC )
  -d, --direction=...     use dipole orientation (001, 111, or 011)
  -c, --constant=...      use the specified lattice constant
  -n                      use the specified number of unit cells
  -o, --output-file=...   use specified output (.xyz) file

Type "A" arrays have nearest neighbor strings of antiparallel dipoles.

Type "B" arrays have nearest neighbor strings of antiparallel dipoles
if the dipoles are contained in a plane perpendicular to the dipole
direction that passes through the dipole.

Example:
   buildDipolarArray -a -l fcc -d 001 -c 5 -n 3 -o A_fcc_001.xyz

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Rev$"
__date__ = "$LastChangedDate$"

__copyright__ = "Copyright (c) 2013 by the University of Notre Dame"
__license__ = "OpenMD"

import sys
import getopt
import string
import math
import numpy

def usage():
    print __doc__

    
    
def createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName):
    # The following section creates 24 basic arrays from Luttinger and
    # Tisza:

    # The six unit vectors are: 3 spatial and 3 to describe the
    # orientation of the dipole.
    
    e1 = numpy.array([1.0,0.0,0.0])
    e2 = numpy.array([0.0,1.0,0.0])
    e3 = numpy.array([0.0,0.0,1.0])

    # Parameters describing the 8 basic arrays:
    cell = numpy.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],
                        [1,1,0],[0,1,1],[1,0,1],[1,1,1]])
    # order in which the basic arrays are constructed in the l loops below:
    corners = numpy.array([[0,0,0],[0,0,1],[0,1,0],[0,1,1],
                           [1,0,0],[1,0,1],[1,1,0],[1,1,1]])

    X = numpy.zeros(192).reshape((8,8,3))
    Y = numpy.zeros(192).reshape((8,8,3))
    Z = numpy.zeros(192).reshape((8,8,3))

    # create the 24 basic arrays using Eq. 12 in Luttinger & Tisza:
    for i in range(8):
        which = 0
        for l1 in range(2):
            for l2 in range(2):
                for l3 in range(2):
                    lvals = numpy.array([l1,l2,l3])
                    value = math.pow(-1, numpy.dot(cell[i], lvals))
                    X[i][which] = value * e1
                    Y[i][which] = value * e2
                    Z[i][which] = value * e3
                    which = which+1

    
    lp_array = numpy.zeros(0).reshape((0,3))
    for i in range(8):        
        lp_array = numpy.vstack((lp_array, corners[i]))
    
    bc_array = numpy.zeros(0).reshape((0,3))
    for i in range(8):
        bc_array = numpy.vstack((bc_array, corners[i] + [0.5,0.5,0.5]))

    xy_array = numpy.zeros(0).reshape((0,3))
    for i in range(8):
        xy_array = numpy.vstack((xy_array, corners[i] + [0.5,0.5,0.0]))

    xz_array = numpy.zeros(0).reshape((0,3))
    for i in range(8):
        xz_array = numpy.vstack((xz_array, corners[i] + [0.5,0.0,0.5]))

    yz_array = numpy.zeros(0).reshape((0,3))
    for i in range(8):
        yz_array = numpy.vstack((yz_array, corners[i] + [0.0,0.5,0.5]))

    known_case = False
    basic_array = numpy.zeros(0).reshape((0,3,3))

    lp_part = numpy.zeros(0).reshape((0,3,3))
    bc_part = numpy.zeros(0).reshape((0,3,3))
    xy_part = numpy.zeros(0).reshape((0,3,3))
    xz_part = numpy.zeros(0).reshape((0,3,3))
    yz_part = numpy.zeros(0).reshape((0,3,3))

    if (arrayType == 'X'):
        if (int(latticeType)):
            which = int(latticeType) - 1
            basic_array = numpy.append(lp_array, X[which], axis=1)
            known_case = True
    if (arrayType == 'Y'):
        if (int(latticeType)):
            which = int(latticeType) - 1
            basic_array = numpy.append(lp_array, Y[which], axis=1)
            known_case = True
    if (arrayType == 'Z'):
        if (int(latticeType)):
            which = int(latticeType) - 1 
            basic_array = numpy.append(lp_array, Z[which], axis=1)
            known_case = True
    if (arrayType == 'A'):
        if (latticeType.lower() == 'sc'):
            basic_array = numpy.append(lp_array, Z[4], axis=1)
            known_case = True
        if (latticeType.lower() == 'bcc'):
            if (dipoleDirection.lower() == '001'):
                lp_part = numpy.append(lp_array,  Z[0], axis=1)
                bc_part = numpy.append(bc_array, -Z[0], axis=1)
                basic_array = numpy.append(lp_part, bc_part, axis=0)
                known_case = True
            if (dipoleDirection.lower() == '111'):
                lp_part = numpy.append(lp_array, Z[4]+X[6]+Y[5], axis=1)
                bc_part = numpy.append(bc_array, Z[4]+X[6]+Y[5], axis=1)
                basic_array = numpy.append(lp_part, bc_part, axis=0)
                known_case = True
        if (latticeType.lower() == 'fcc'):
            if (dipoleDirection.lower() == '001'):
                lp_part = numpy.append(lp_array,  Z[0], axis=1)
                xy_part = numpy.append(xy_array,  Z[0], axis=1)
                xz_part = numpy.append(xz_array, -Z[0], axis=1)
                yz_part = numpy.append(yz_array, -Z[0], axis=1)

                basic_array = numpy.append(lp_part, xy_part, axis=0)
                basic_array = numpy.append(basic_array, xz_part, axis=0)
                basic_array = numpy.append(basic_array, yz_part, axis=0)
                
                known_case = True
            if (dipoleDirection.lower() == '011'):
                lp_part = numpy.append(lp_array,  Z[0]+Y[0], axis=1)
                xy_part = numpy.append(xy_array, -Z[0]-Y[0], axis=1)
                xz_part = numpy.append(xz_array, -Z[0]-Y[0], axis=1)
                yz_part = numpy.append(yz_array,  Z[0]+Y[0], axis=1)

                basic_array = numpy.append(lp_part, xy_part, axis=0)
                basic_array = numpy.append(basic_array, xz_part, axis=0)
                basic_array = numpy.append(basic_array, yz_part, axis=0)

                known_case = True
    else:
        if (latticeType.lower() == 'sc'):
            basic_array = numpy.append(lp_array, Z[4], axis=1)
            known_case = True
        if (latticeType.lower() == 'bcc'):
            if (dipoleDirection.lower() == '001'):
                lp_part = numpy.append(lp_array,  Z[4], axis=1)
                bc_part = numpy.append(bc_array, -Z[4], axis=1)
                basic_array = numpy.append(lp_part, bc_part, axis=0)
                known_case = True
            if (dipoleDirection.lower() == '111'):
                lp_part = numpy.append(lp_array, Z[4]+X[6]+Y[5], axis=1)
                bc_part = numpy.append(bc_array, Z[4]+X[6]+Y[5], axis=1)
                basic_array = numpy.append(lp_part, bc_part, axis=0)
                known_case = True
        if (latticeType.lower() == 'fcc'):
            if (dipoleDirection.lower() == '001'):
                lp_part = numpy.append(lp_array,  Z[0], axis=1)
                xy_part = numpy.append(xy_array, -Z[0], axis=1)
                xz_part = numpy.append(xz_array, -Z[0], axis=1)
                yz_part = numpy.append(yz_array,  Z[0], axis=1)

                basic_array = numpy.append(lp_part, xy_part, axis=0)
                basic_array = numpy.append(basic_array, xz_part, axis=0)
                basic_array = numpy.append(basic_array, yz_part, axis=0)

                known_case = True
            if (dipoleDirection.lower() == '011'):
                lp_part = numpy.append(lp_array,  Z[7]+Y[7], axis=1)
                xy_part = numpy.append(xy_array,  Z[7]+Y[7], axis=1)
                xz_part = numpy.append(xz_array, -Z[7]-Y[7], axis=1)
                yz_part = numpy.append(yz_array,  Z[7]+Y[7], axis=1)

                basic_array = numpy.append(lp_part, xy_part, axis=0)
                basic_array = numpy.append(basic_array, xz_part, axis=0)
                basic_array = numpy.append(basic_array, yz_part, axis=0)

                known_case = True

    if (not known_case):
        print "unhandled combination of lattice and dipole direction"
        print __doc__

    bravais_lattice = numpy.zeros(0).reshape((0,6))
    for i in range(latticeNumber):
        for j in range(latticeNumber):
            for k in range(latticeNumber):
                for l in range(len(basic_array)):
                    lat_vec = numpy.array([[2*i, 2*j, 2*k, 0.0, 0.0, 0.0]])
                    bravais_lattice = numpy.append(bravais_lattice, lat_vec + basic_array[l], axis=0)

    outputFile = open(outputFileName, 'w')
    
    outputFile.write('<OpenMD version=2>\n')
    outputFile.write('  <MetaData>\n')
    outputFile.write('     molecule{\n')
    outputFile.write('       name = \"D\";\n')
    outputFile.write('       \n')
    outputFile.write('       atom[0]{\n')
    outputFile.write('         type = \"D\";\n')
    outputFile.write('         position(0.0, 0.0, 0.0);\n')
    outputFile.write('         orientation(0.0, 0.0, 0.0);\n')
    outputFile.write('       }\n')
    outputFile.write('     }\n')
    outputFile.write('     component{\n')
    outputFile.write('       type = \"D\";\n')
    outputFile.write('       nMol = '+ repr(len(bravais_lattice)) + ';\n')
    outputFile.write('     }\n')

    outputFile.write('     ensemble = NVE;\n')
    outputFile.write('     forceField = \"Multipole\";\n')

    outputFile.write('     cutoffMethod = \"shifted_force\";\n')
    outputFile.write('     electrostaticScreeningMethod = \"damped\";\n')

    outputFile.write('     cutoffRadius = 9.0;\n')
    outputFile.write('     dampingAlpha = 0.18;\n')
    outputFile.write('     statFileFormat = \"TIME|TOTAL_ENERGY|POTENTIAL_ENERGY|KINETIC_ENERGY|TEMPERATURE|PRESSURE|VOLUME|CONSERVED_QUANTITY|ELECTROSTATIC_POTENTIAL\";\n')
    outputFile.write('     dt = 1.0;\n')
    outputFile.write('     runTime = 1.0;\n')
    outputFile.write('     sampleTime = 1.0;\n')
    outputFile.write('     statusTime = 1.0;\n')
    outputFile.write('  </MetaData>\n')
    outputFile.write('  <Snapshot>\n')
    outputFile.write('    <FrameData>\n');
    outputFile.write("        Time: %.10g\n" % (0.0))
    
    Hxx = 2.0 * latticeConstant * latticeNumber
    Hyy = 2.0 * latticeConstant * latticeNumber
    Hzz = 2.0 * latticeConstant * latticeNumber
    
    outputFile.write('        Hmat: {{%d, 0, 0}, {0, %d, 0}, {0, 0, %d}}\n' % (Hxx, Hyy, Hzz))
    outputFile.write('    </FrameData>\n')
    outputFile.write('    <StuntDoubles>\n')
    sdFormat = 'pvqj'
    index = 0

    for i in range(len(bravais_lattice)):
        xcart = latticeConstant*(bravais_lattice[i][0])
        ycart = latticeConstant*(bravais_lattice[i][1])
        zcart = latticeConstant*(bravais_lattice[i][2])
        dx = bravais_lattice[i][3]
        dy = bravais_lattice[i][4]
        dz = bravais_lattice[i][5]

        dlen = math.sqrt(dx*dx + dy*dy + dz*dz)
        ctheta = dz / dlen
        theta = math.acos(ctheta)
        stheta = math.sqrt(1.0 - ctheta*ctheta)
        psi = 0.0
        phi = math.atan2(-dy/dlen, dx/dlen)

        q = [0.0,0.0,0.0,0.0]
        q[0] = math.cos(theta/2)*math.cos((phi+psi)/2)
        q[1] = math.sin(theta/2)*math.cos((phi-psi)/2)
        q[2] = math.sin(theta/2)*math.sin((phi-psi)/2)
        q[3] = math.cos(theta/2)*math.sin((phi+psi)/2)

        qlen = math.sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3])
        q[0] = q[0]/qlen
        q[1] = q[1]/qlen
        q[2] = q[2]/qlen
        q[3] = q[3]/qlen
        
        outputFile.write("%10d %7s %g %g %1g %g %g %g %13e %13e %13e %13e %g %g %g\n" % (index, sdFormat, xcart, ycart, zcart, 0.0, 0.0, 0.0, q[0], q[1], q[2], q[3], 0.0, 0.0, 0.0))
        index = index+1

    outputFile.write("    </StuntDoubles>\n")
    outputFile.write("  </Snapshot>\n")
    outputFile.write("</OpenMD>\n")
    outputFile.close()

    outputFile.close()
    
def main(argv):
    
    arrayType = "A"
    haveOutputFileName = False
    latticeType = "fcc"
    dipoleDirection = "001"
    latticeNumber = 3
    latticeConstant = 4
    try:                                
        opts, args = getopt.getopt(argv, "hxyzabl:d:c:n:o:", ["help","array-type-X", "array-type-Y", "array-type-Z", "array-type-A", "array-type-B", "lattice=" "direction=", "constant=", "output-file="])  
    except getopt.GetoptError:           
        usage()                          
        sys.exit(2)                     
    for opt, arg in opts:                
        if opt in ("-h", "--help"):      
            usage()                     
            sys.exit()
        elif opt in ("-x", "--array-type-X"):
            arrayType = "X"
        elif opt in ("-y", "--array-type-Y"):
            arrayType = "Y"
        elif opt in ("-z", "--array-type-Z"):
            arrayType = "Z"
        elif opt in ("-b", "--array-type-B"):
            arrayType = "B"
        elif opt in ("-l", "--lattice"): 
            latticeType = arg
        elif opt in ("-d", "--direction"):
            dipoleDirection = arg
        elif opt in ("-c", "--constant"):
            latticeConstant = float(arg)
        elif opt in ("-n"):
            latticeNumber = int(arg)
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            haveOutputFileName = True
    if (not haveOutputFileName):
        usage()
        print "No output file was specified"
        sys.exit()
        
    createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName);

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1911
Committed:	Tue Jul 23 01:58:19 2013 UTC (12 years, 1 month ago) by gezelter
File size:	13965 byte(s)
Log Message:	modified svn propset info
#	User	Rev	Content
1	gezelter	1889	#! /usr/bin/env python
2
3			"""Dipolar Lattice Builder
4
5			Creates cubic lattices of dipoles to test the dipole-dipole
6			interaction code.
7
8			Usage: buildDipolarArray
9
10			Options:
11			-h, --help show this help
12	gezelter	1899	-x, --array-type-X use one of the basic "X" arrays
13			-y, --array-type-Y use one of the basic "Y" arrays
14			-z, --array-type-Z use one of the basic "Z" arrays
15	gezelter	1889	-a, --array-type-A use array type "A" (default)
16			-b, --array-type-B use array type "B"
17			-l, --lattice=... use the specified lattice ( SC, FCC, or BCC )
18			-d, --direction=... use dipole orientation (001, 111, or 011)
19			-c, --constant=... use the specified lattice constant
20			-n use the specified number of unit cells
21			-o, --output-file=... use specified output (.xyz) file
22
23			Type "A" arrays have nearest neighbor strings of antiparallel dipoles.
24
25			Type "B" arrays have nearest neighbor strings of antiparallel dipoles
26			if the dipoles are contained in a plane perpendicular to the dipole
27			direction that passes through the dipole.
28
29			Example:
30			buildDipolarArray -a -l fcc -d 001 -c 5 -n 3 -o A_fcc_001.xyz
31
32			"""
33
34			__author__ = "Dan Gezelter (gezelter@nd.edu)"
35	gezelter	1911	__version__ = "$Rev$"
36			__date__ = "$LastChangedDate$"
37	gezelter	1889
38			__copyright__ = "Copyright (c) 2013 by the University of Notre Dame"
39			__license__ = "OpenMD"
40
41			import sys
42			import getopt
43			import string
44			import math
45			import numpy
46
47			def usage():
48			print __doc__
49	gezelter	1898
50	gezelter	1889
51	gezelter	1898
52	gezelter	1899	def createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName):
53	gezelter	1898	# The following section creates 24 basic arrays from Luttinger and
54			# Tisza:
55	gezelter	1889
56	gezelter	1898	# The six unit vectors are: 3 spatial and 3 to describe the
57			# orientation of the dipole.
58
59	gezelter	1907	e1 = numpy.array([1.0,0.0,0.0])
60			e2 = numpy.array([0.0,1.0,0.0])
61			e3 = numpy.array([0.0,0.0,1.0])
62	gezelter	1898
63			# Parameters describing the 8 basic arrays:
64			cell = numpy.array([[0,0,0],[0,0,1],[1,0,0],[0,1,0],
65			[1,1,0],[0,1,1],[1,0,1],[1,1,1]])
66	gezelter	1909	# order in which the basic arrays are constructed in the l loops below:
67			corners = numpy.array([[0,0,0],[0,0,1],[0,1,0],[0,1,1],
68			[1,0,0],[1,0,1],[1,1,0],[1,1,1]])
69	gezelter	1898
70	gezelter	1907	X = numpy.zeros(192).reshape((8,8,3))
71			Y = numpy.zeros(192).reshape((8,8,3))
72			Z = numpy.zeros(192).reshape((8,8,3))
73	gezelter	1898
74			# create the 24 basic arrays using Eq. 12 in Luttinger & Tisza:
75			for i in range(8):
76			which = 0
77			for l1 in range(2):
78			for l2 in range(2):
79			for l3 in range(2):
80			lvals = numpy.array([l1,l2,l3])
81			value = math.pow(-1, numpy.dot(cell[i], lvals))
82	gezelter	1907	X[i][which] = value * e1
83			Y[i][which] = value * e2
84			Z[i][which] = value * e3
85	gezelter	1909	which = which+1
86
87	gezelter	1907
88			lp_array = numpy.zeros(0).reshape((0,3))
89			for i in range(8):
90	gezelter	1909	lp_array = numpy.vstack((lp_array, corners[i]))
91	gezelter	1907
92			bc_array = numpy.zeros(0).reshape((0,3))
93			for i in range(8):
94	gezelter	1909	bc_array = numpy.vstack((bc_array, corners[i] + [0.5,0.5,0.5]))
95	gezelter	1898
96	gezelter	1907	xy_array = numpy.zeros(0).reshape((0,3))
97			for i in range(8):
98	gezelter	1909	xy_array = numpy.vstack((xy_array, corners[i] + [0.5,0.5,0.0]))
99	gezelter	1898
100	gezelter	1907	xz_array = numpy.zeros(0).reshape((0,3))
101			for i in range(8):
102	gezelter	1909	xz_array = numpy.vstack((xz_array, corners[i] + [0.5,0.0,0.5]))
103	gezelter	1898
104	gezelter	1907	yz_array = numpy.zeros(0).reshape((0,3))
105			for i in range(8):
106	gezelter	1909	yz_array = numpy.vstack((yz_array, corners[i] + [0.0,0.5,0.5]))
107	gezelter	1898
108	gezelter	1907	known_case = False
109			basic_array = numpy.zeros(0).reshape((0,3,3))
110	gezelter	1889
111	gezelter	1907	lp_part = numpy.zeros(0).reshape((0,3,3))
112			bc_part = numpy.zeros(0).reshape((0,3,3))
113			xy_part = numpy.zeros(0).reshape((0,3,3))
114			xz_part = numpy.zeros(0).reshape((0,3,3))
115			yz_part = numpy.zeros(0).reshape((0,3,3))
116	gezelter	1889
117	gezelter	1899	if (arrayType == 'X'):
118			if (int(latticeType)):
119			which = int(latticeType) - 1
120	gezelter	1907	basic_array = numpy.append(lp_array, X[which], axis=1)
121	gezelter	1899	known_case = True
122			if (arrayType == 'Y'):
123			if (int(latticeType)):
124			which = int(latticeType) - 1
125	gezelter	1907	basic_array = numpy.append(lp_array, Y[which], axis=1)
126	gezelter	1899	known_case = True
127			if (arrayType == 'Z'):
128			if (int(latticeType)):
129			which = int(latticeType) - 1
130	gezelter	1907	basic_array = numpy.append(lp_array, Z[which], axis=1)
131	gezelter	1899	known_case = True
132			if (arrayType == 'A'):
133	gezelter	1898	if (latticeType.lower() == 'sc'):
134	gezelter	1907	basic_array = numpy.append(lp_array, Z[4], axis=1)
135	gezelter	1898	known_case = True
136			if (latticeType.lower() == 'bcc'):
137			if (dipoleDirection.lower() == '001'):
138	gezelter	1907	lp_part = numpy.append(lp_array, Z[0], axis=1)
139			bc_part = numpy.append(bc_array, -Z[0], axis=1)
140			basic_array = numpy.append(lp_part, bc_part, axis=0)
141	gezelter	1898	known_case = True
142			if (dipoleDirection.lower() == '111'):
143	gezelter	1907	lp_part = numpy.append(lp_array, Z[4]+X[6]+Y[5], axis=1)
144			bc_part = numpy.append(bc_array, Z[4]+X[6]+Y[5], axis=1)
145			basic_array = numpy.append(lp_part, bc_part, axis=0)
146	gezelter	1898	known_case = True
147			if (latticeType.lower() == 'fcc'):
148			if (dipoleDirection.lower() == '001'):
149	gezelter	1907	lp_part = numpy.append(lp_array, Z[0], axis=1)
150			xy_part = numpy.append(xy_array, Z[0], axis=1)
151			xz_part = numpy.append(xz_array, -Z[0], axis=1)
152			yz_part = numpy.append(yz_array, -Z[0], axis=1)
153
154			basic_array = numpy.append(lp_part, xy_part, axis=0)
155			basic_array = numpy.append(basic_array, xz_part, axis=0)
156			basic_array = numpy.append(basic_array, yz_part, axis=0)
157
158	gezelter	1898	known_case = True
159			if (dipoleDirection.lower() == '011'):
160	gezelter	1907	lp_part = numpy.append(lp_array, Z[0]+Y[0], axis=1)
161			xy_part = numpy.append(xy_array, -Z[0]-Y[0], axis=1)
162			xz_part = numpy.append(xz_array, -Z[0]-Y[0], axis=1)
163			yz_part = numpy.append(yz_array, Z[0]+Y[0], axis=1)
164
165			basic_array = numpy.append(lp_part, xy_part, axis=0)
166			basic_array = numpy.append(basic_array, xz_part, axis=0)
167			basic_array = numpy.append(basic_array, yz_part, axis=0)
168
169	gezelter	1898	known_case = True
170	gezelter	1889	else:
171			if (latticeType.lower() == 'sc'):
172	gezelter	1907	basic_array = numpy.append(lp_array, Z[4], axis=1)
173	gezelter	1898	known_case = True
174			if (latticeType.lower() == 'bcc'):
175	gezelter	1889	if (dipoleDirection.lower() == '001'):
176	gezelter	1907	lp_part = numpy.append(lp_array, Z[4], axis=1)
177			bc_part = numpy.append(bc_array, -Z[4], axis=1)
178			basic_array = numpy.append(lp_part, bc_part, axis=0)
179	gezelter	1898	known_case = True
180			if (dipoleDirection.lower() == '111'):
181	gezelter	1907	lp_part = numpy.append(lp_array, Z[4]+X[6]+Y[5], axis=1)
182			bc_part = numpy.append(bc_array, Z[4]+X[6]+Y[5], axis=1)
183			basic_array = numpy.append(lp_part, bc_part, axis=0)
184	gezelter	1898	known_case = True
185			if (latticeType.lower() == 'fcc'):
186	gezelter	1889	if (dipoleDirection.lower() == '001'):
187	gezelter	1907	lp_part = numpy.append(lp_array, Z[0], axis=1)
188			xy_part = numpy.append(xy_array, -Z[0], axis=1)
189			xz_part = numpy.append(xz_array, -Z[0], axis=1)
190			yz_part = numpy.append(yz_array, Z[0], axis=1)
191
192			basic_array = numpy.append(lp_part, xy_part, axis=0)
193			basic_array = numpy.append(basic_array, xz_part, axis=0)
194			basic_array = numpy.append(basic_array, yz_part, axis=0)
195
196	gezelter	1898	known_case = True
197			if (dipoleDirection.lower() == '011'):
198	gezelter	1907	lp_part = numpy.append(lp_array, Z[7]+Y[7], axis=1)
199			xy_part = numpy.append(xy_array, Z[7]+Y[7], axis=1)
200			xz_part = numpy.append(xz_array, -Z[7]-Y[7], axis=1)
201			yz_part = numpy.append(yz_array, Z[7]+Y[7], axis=1)
202
203			basic_array = numpy.append(lp_part, xy_part, axis=0)
204			basic_array = numpy.append(basic_array, xz_part, axis=0)
205			basic_array = numpy.append(basic_array, yz_part, axis=0)
206
207	gezelter	1898	known_case = True
208	gezelter	1889
209	gezelter	1898	if (not known_case):
210			print "unhandled combination of lattice and dipole direction"
211			print __doc__
212
213	gezelter	1907	bravais_lattice = numpy.zeros(0).reshape((0,6))
214	gezelter	1889	for i in range(latticeNumber):
215			for j in range(latticeNumber):
216			for k in range(latticeNumber):
217	gezelter	1898	for l in range(len(basic_array)):
218	gezelter	1907	lat_vec = numpy.array([[2i, 2j, 2*k, 0.0, 0.0, 0.0]])
219			bravais_lattice = numpy.append(bravais_lattice, lat_vec + basic_array[l], axis=0)
220	gezelter	1889
221	gezelter	1898	outputFile = open(outputFileName, 'w')
222
223			outputFile.write('<OpenMD version=2>\n')
224			outputFile.write(' <MetaData>\n')
225			outputFile.write(' molecule{\n')
226			outputFile.write(' name = \"D\";\n')
227			outputFile.write(' \n')
228			outputFile.write(' atom[0]{\n')
229			outputFile.write(' type = \"D\";\n')
230			outputFile.write(' position(0.0, 0.0, 0.0);\n')
231			outputFile.write(' orientation(0.0, 0.0, 0.0);\n')
232			outputFile.write(' }\n')
233			outputFile.write(' }\n')
234			outputFile.write(' component{\n')
235			outputFile.write(' type = \"D\";\n')
236			outputFile.write(' nMol = '+ repr(len(bravais_lattice)) + ';\n')
237			outputFile.write(' }\n')
238	gezelter	1889
239	gezelter	1898	outputFile.write(' ensemble = NVE;\n')
240			outputFile.write(' forceField = \"Multipole\";\n')
241
242			outputFile.write(' cutoffMethod = \"shifted_force\";\n')
243			outputFile.write(' electrostaticScreeningMethod = \"damped\";\n')
244
245			outputFile.write(' cutoffRadius = 9.0;\n')
246			outputFile.write(' dampingAlpha = 0.18;\n')
247			outputFile.write(' statFileFormat = \"TIME\|TOTAL_ENERGY\|POTENTIAL_ENERGY\|KINETIC_ENERGY\|TEMPERATURE\|PRESSURE\|VOLUME\|CONSERVED_QUANTITY\|ELECTROSTATIC_POTENTIAL\";\n')
248			outputFile.write(' dt = 1.0;\n')
249			outputFile.write(' runTime = 1.0;\n')
250			outputFile.write(' sampleTime = 1.0;\n')
251			outputFile.write(' statusTime = 1.0;\n')
252			outputFile.write(' </MetaData>\n')
253			outputFile.write(' <Snapshot>\n')
254			outputFile.write(' <FrameData>\n');
255			outputFile.write(" Time: %.10g\n" % (0.0))
256	gezelter	1889
257	gezelter	1898	Hxx = 2.0 * latticeConstant * latticeNumber
258			Hyy = 2.0 * latticeConstant * latticeNumber
259			Hzz = 2.0 * latticeConstant * latticeNumber
260
261			outputFile.write(' Hmat: {{%d, 0, 0}, {0, %d, 0}, {0, 0, %d}}\n' % (Hxx, Hyy, Hzz))
262			outputFile.write(' </FrameData>\n')
263			outputFile.write(' <StuntDoubles>\n')
264			sdFormat = 'pvqj'
265			index = 0
266	gezelter	1907
267	gezelter	1889	for i in range(len(bravais_lattice)):
268			xcart = latticeConstant*(bravais_lattice[i][0])
269			ycart = latticeConstant*(bravais_lattice[i][1])
270			zcart = latticeConstant*(bravais_lattice[i][2])
271	gezelter	1898	dx = bravais_lattice[i][3]
272			dy = bravais_lattice[i][4]
273			dz = bravais_lattice[i][5]
274
275	gezelter	1909	dlen = math.sqrt(dxdx + dydy + dz*dz)
276			ctheta = dz / dlen
277			theta = math.acos(ctheta)
278			stheta = math.sqrt(1.0 - ctheta*ctheta)
279			psi = 0.0
280			phi = math.atan2(-dy/dlen, dx/dlen)
281	gezelter	1898
282	gezelter	1909	q = [0.0,0.0,0.0,0.0]
283			q[0] = math.cos(theta/2)*math.cos((phi+psi)/2)
284			q[1] = math.sin(theta/2)*math.cos((phi-psi)/2)
285			q[2] = math.sin(theta/2)*math.sin((phi-psi)/2)
286			q[3] = math.cos(theta/2)*math.sin((phi+psi)/2)
287	gezelter	1898
288	gezelter	1905	qlen = math.sqrt(q[0]q[0] + q[1]q[1] + q[2]q[2] + q[3]q[3])
289			q[0] = q[0]/qlen
290			q[1] = q[1]/qlen
291			q[2] = q[2]/qlen
292			q[3] = q[3]/qlen
293	gezelter	1898
294	gezelter	1909	outputFile.write("%10d %7s %g %g %1g %g %g %g %13e %13e %13e %13e %g %g %g\n" % (index, sdFormat, xcart, ycart, zcart, 0.0, 0.0, 0.0, q[0], q[1], q[2], q[3], 0.0, 0.0, 0.0))
295	gezelter	1898	index = index+1
296
297			outputFile.write(" </StuntDoubles>\n")
298			outputFile.write(" </Snapshot>\n")
299			outputFile.write("</OpenMD>\n")
300	gezelter	1889	outputFile.close()
301
302	gezelter	1898	outputFile.close()
303
304	gezelter	1889	def main(argv):
305
306	gezelter	1899	arrayType = "A"
307	gezelter	1889	haveOutputFileName = False
308			latticeType = "fcc"
309			dipoleDirection = "001"
310			latticeNumber = 3
311			latticeConstant = 4
312			try:
313	gezelter	1899	opts, args = getopt.getopt(argv, "hxyzabl:d:c:n:o:", ["help","array-type-X", "array-type-Y", "array-type-Z", "array-type-A", "array-type-B", "lattice=" "direction=", "constant=", "output-file="])
314	gezelter	1889	except getopt.GetoptError:
315			usage()
316			sys.exit(2)
317			for opt, arg in opts:
318			if opt in ("-h", "--help"):
319			usage()
320			sys.exit()
321	gezelter	1899	elif opt in ("-x", "--array-type-X"):
322			arrayType = "X"
323			elif opt in ("-y", "--array-type-Y"):
324			arrayType = "Y"
325			elif opt in ("-z", "--array-type-Z"):
326			arrayType = "Z"
327	gezelter	1889	elif opt in ("-b", "--array-type-B"):
328	gezelter	1899	arrayType = "B"
329	gezelter	1889	elif opt in ("-l", "--lattice"):
330			latticeType = arg
331			elif opt in ("-d", "--direction"):
332			dipoleDirection = arg
333			elif opt in ("-c", "--constant"):
334			latticeConstant = float(arg)
335			elif opt in ("-n"):
336			latticeNumber = int(arg)
337			elif opt in ("-o", "--output-file"):
338			outputFileName = arg
339			haveOutputFileName = True
340			if (not haveOutputFileName):
341			usage()
342			print "No output file was specified"
343			sys.exit()
344
345	gezelter	1899	createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName);
346	gezelter	1889
347			if __name__ == "__main__":
348			if len(sys.argv) == 1:
349			usage()
350			sys.exit()
351			main(sys.argv[1:])