[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 (dipoleDirection.lower() == 'min'):
                lp_part = numpy.append(lp_array, X[6]+Y[4]+Z[5], axis=1)
                bc_part = numpy.append(bc_array, X[4]+Y[5]+Z[6], 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(dx/dlen, -dy/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:	1914
Committed:	Mon Jul 29 15:34:04 2013 UTC (12 years, 1 month ago) by gezelter
File size:	14265 byte(s)
Log Message:	Fixed a quaternion bug, updated the affected lattices.
#	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	gezelter	1912	if (dipoleDirection.lower() == 'min'):
148	gezelter	1914	lp_part = numpy.append(lp_array, X[6]+Y[4]+Z[5], axis=1)
149			bc_part = numpy.append(bc_array, X[4]+Y[5]+Z[6], axis=1)
150	gezelter	1912	basic_array = numpy.append(lp_part, bc_part, axis=0)
151			known_case = True
152	gezelter	1898	if (latticeType.lower() == 'fcc'):
153			if (dipoleDirection.lower() == '001'):
154	gezelter	1907	lp_part = numpy.append(lp_array, Z[0], axis=1)
155			xy_part = numpy.append(xy_array, Z[0], axis=1)
156			xz_part = numpy.append(xz_array, -Z[0], axis=1)
157			yz_part = numpy.append(yz_array, -Z[0], axis=1)
158
159			basic_array = numpy.append(lp_part, xy_part, axis=0)
160			basic_array = numpy.append(basic_array, xz_part, axis=0)
161			basic_array = numpy.append(basic_array, yz_part, axis=0)
162
163	gezelter	1898	known_case = True
164			if (dipoleDirection.lower() == '011'):
165	gezelter	1907	lp_part = numpy.append(lp_array, Z[0]+Y[0], axis=1)
166			xy_part = numpy.append(xy_array, -Z[0]-Y[0], axis=1)
167			xz_part = numpy.append(xz_array, -Z[0]-Y[0], axis=1)
168			yz_part = numpy.append(yz_array, Z[0]+Y[0], axis=1)
169
170			basic_array = numpy.append(lp_part, xy_part, axis=0)
171			basic_array = numpy.append(basic_array, xz_part, axis=0)
172			basic_array = numpy.append(basic_array, yz_part, axis=0)
173
174	gezelter	1898	known_case = True
175	gezelter	1889	else:
176			if (latticeType.lower() == 'sc'):
177	gezelter	1907	basic_array = numpy.append(lp_array, Z[4], axis=1)
178	gezelter	1898	known_case = True
179			if (latticeType.lower() == 'bcc'):
180	gezelter	1889	if (dipoleDirection.lower() == '001'):
181	gezelter	1907	lp_part = numpy.append(lp_array, Z[4], axis=1)
182			bc_part = numpy.append(bc_array, -Z[4], axis=1)
183			basic_array = numpy.append(lp_part, bc_part, axis=0)
184	gezelter	1898	known_case = True
185			if (dipoleDirection.lower() == '111'):
186	gezelter	1907	lp_part = numpy.append(lp_array, Z[4]+X[6]+Y[5], axis=1)
187			bc_part = numpy.append(bc_array, Z[4]+X[6]+Y[5], axis=1)
188			basic_array = numpy.append(lp_part, bc_part, axis=0)
189	gezelter	1898	known_case = True
190			if (latticeType.lower() == 'fcc'):
191	gezelter	1889	if (dipoleDirection.lower() == '001'):
192	gezelter	1907	lp_part = numpy.append(lp_array, Z[0], axis=1)
193			xy_part = numpy.append(xy_array, -Z[0], axis=1)
194			xz_part = numpy.append(xz_array, -Z[0], axis=1)
195			yz_part = numpy.append(yz_array, Z[0], axis=1)
196
197			basic_array = numpy.append(lp_part, xy_part, axis=0)
198			basic_array = numpy.append(basic_array, xz_part, axis=0)
199			basic_array = numpy.append(basic_array, yz_part, axis=0)
200
201	gezelter	1898	known_case = True
202			if (dipoleDirection.lower() == '011'):
203	gezelter	1907	lp_part = numpy.append(lp_array, Z[7]+Y[7], axis=1)
204			xy_part = numpy.append(xy_array, Z[7]+Y[7], axis=1)
205			xz_part = numpy.append(xz_array, -Z[7]-Y[7], axis=1)
206			yz_part = numpy.append(yz_array, Z[7]+Y[7], axis=1)
207
208			basic_array = numpy.append(lp_part, xy_part, axis=0)
209			basic_array = numpy.append(basic_array, xz_part, axis=0)
210			basic_array = numpy.append(basic_array, yz_part, axis=0)
211
212	gezelter	1898	known_case = True
213	gezelter	1889
214	gezelter	1898	if (not known_case):
215			print "unhandled combination of lattice and dipole direction"
216			print __doc__
217
218	gezelter	1907	bravais_lattice = numpy.zeros(0).reshape((0,6))
219	gezelter	1889	for i in range(latticeNumber):
220			for j in range(latticeNumber):
221			for k in range(latticeNumber):
222	gezelter	1898	for l in range(len(basic_array)):
223	gezelter	1907	lat_vec = numpy.array([[2i, 2j, 2*k, 0.0, 0.0, 0.0]])
224			bravais_lattice = numpy.append(bravais_lattice, lat_vec + basic_array[l], axis=0)
225	gezelter	1889
226	gezelter	1898	outputFile = open(outputFileName, 'w')
227
228			outputFile.write('<OpenMD version=2>\n')
229			outputFile.write(' <MetaData>\n')
230			outputFile.write(' molecule{\n')
231			outputFile.write(' name = \"D\";\n')
232			outputFile.write(' \n')
233			outputFile.write(' atom[0]{\n')
234			outputFile.write(' type = \"D\";\n')
235			outputFile.write(' position(0.0, 0.0, 0.0);\n')
236			outputFile.write(' orientation(0.0, 0.0, 0.0);\n')
237			outputFile.write(' }\n')
238			outputFile.write(' }\n')
239			outputFile.write(' component{\n')
240			outputFile.write(' type = \"D\";\n')
241			outputFile.write(' nMol = '+ repr(len(bravais_lattice)) + ';\n')
242			outputFile.write(' }\n')
243	gezelter	1889
244	gezelter	1898	outputFile.write(' ensemble = NVE;\n')
245			outputFile.write(' forceField = \"Multipole\";\n')
246
247			outputFile.write(' cutoffMethod = \"shifted_force\";\n')
248			outputFile.write(' electrostaticScreeningMethod = \"damped\";\n')
249
250			outputFile.write(' cutoffRadius = 9.0;\n')
251			outputFile.write(' dampingAlpha = 0.18;\n')
252			outputFile.write(' statFileFormat = \"TIME\|TOTAL_ENERGY\|POTENTIAL_ENERGY\|KINETIC_ENERGY\|TEMPERATURE\|PRESSURE\|VOLUME\|CONSERVED_QUANTITY\|ELECTROSTATIC_POTENTIAL\";\n')
253			outputFile.write(' dt = 1.0;\n')
254			outputFile.write(' runTime = 1.0;\n')
255			outputFile.write(' sampleTime = 1.0;\n')
256			outputFile.write(' statusTime = 1.0;\n')
257			outputFile.write(' </MetaData>\n')
258			outputFile.write(' <Snapshot>\n')
259			outputFile.write(' <FrameData>\n');
260			outputFile.write(" Time: %.10g\n" % (0.0))
261	gezelter	1889
262	gezelter	1898	Hxx = 2.0 * latticeConstant * latticeNumber
263			Hyy = 2.0 * latticeConstant * latticeNumber
264			Hzz = 2.0 * latticeConstant * latticeNumber
265
266			outputFile.write(' Hmat: {{%d, 0, 0}, {0, %d, 0}, {0, 0, %d}}\n' % (Hxx, Hyy, Hzz))
267			outputFile.write(' </FrameData>\n')
268			outputFile.write(' <StuntDoubles>\n')
269			sdFormat = 'pvqj'
270			index = 0
271	gezelter	1907
272	gezelter	1889	for i in range(len(bravais_lattice)):
273			xcart = latticeConstant*(bravais_lattice[i][0])
274			ycart = latticeConstant*(bravais_lattice[i][1])
275			zcart = latticeConstant*(bravais_lattice[i][2])
276	gezelter	1898	dx = bravais_lattice[i][3]
277			dy = bravais_lattice[i][4]
278			dz = bravais_lattice[i][5]
279
280	gezelter	1909	dlen = math.sqrt(dxdx + dydy + dz*dz)
281			ctheta = dz / dlen
282			theta = math.acos(ctheta)
283			stheta = math.sqrt(1.0 - ctheta*ctheta)
284			psi = 0.0
285	gezelter	1914	phi = math.atan2(dx/dlen, -dy/dlen)
286	gezelter	1898
287	gezelter	1909	q = [0.0,0.0,0.0,0.0]
288			q[0] = math.cos(theta/2)*math.cos((phi+psi)/2)
289			q[1] = math.sin(theta/2)*math.cos((phi-psi)/2)
290			q[2] = math.sin(theta/2)*math.sin((phi-psi)/2)
291			q[3] = math.cos(theta/2)*math.sin((phi+psi)/2)
292	gezelter	1898
293	gezelter	1905	qlen = math.sqrt(q[0]q[0] + q[1]q[1] + q[2]q[2] + q[3]q[3])
294			q[0] = q[0]/qlen
295			q[1] = q[1]/qlen
296			q[2] = q[2]/qlen
297			q[3] = q[3]/qlen
298	gezelter	1898
299	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))
300	gezelter	1898	index = index+1
301
302			outputFile.write(" </StuntDoubles>\n")
303			outputFile.write(" </Snapshot>\n")
304			outputFile.write("</OpenMD>\n")
305	gezelter	1889	outputFile.close()
306
307	gezelter	1898	outputFile.close()
308
309	gezelter	1889	def main(argv):
310
311	gezelter	1899	arrayType = "A"
312	gezelter	1889	haveOutputFileName = False
313			latticeType = "fcc"
314			dipoleDirection = "001"
315			latticeNumber = 3
316			latticeConstant = 4
317			try:
318	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="])
319	gezelter	1889	except getopt.GetoptError:
320			usage()
321			sys.exit(2)
322			for opt, arg in opts:
323			if opt in ("-h", "--help"):
324			usage()
325			sys.exit()
326	gezelter	1899	elif opt in ("-x", "--array-type-X"):
327			arrayType = "X"
328			elif opt in ("-y", "--array-type-Y"):
329			arrayType = "Y"
330			elif opt in ("-z", "--array-type-Z"):
331			arrayType = "Z"
332	gezelter	1889	elif opt in ("-b", "--array-type-B"):
333	gezelter	1899	arrayType = "B"
334	gezelter	1889	elif opt in ("-l", "--lattice"):
335			latticeType = arg
336			elif opt in ("-d", "--direction"):
337			dipoleDirection = arg
338			elif opt in ("-c", "--constant"):
339			latticeConstant = float(arg)
340			elif opt in ("-n"):
341			latticeNumber = int(arg)
342			elif opt in ("-o", "--output-file"):
343			outputFileName = arg
344			haveOutputFileName = True
345			if (not haveOutputFileName):
346			usage()
347			print "No output file was specified"
348			sys.exit()
349
350	gezelter	1899	createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName);
351	gezelter	1889
352			if __name__ == "__main__":
353			if len(sys.argv) == 1:
354			usage()
355			sys.exit()
356			main(sys.argv[1:])