[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):
        # The X and Y arrays are cubic rotations of the Z array, so we
        # need to re-order them to match the numbering scheme in
        # Luttinger & Tisza:
        if (i > 0 and i < 4):
            iX = 1 + (i + 1) % 3
            iY = 1 + (i ) % 3
        elif (i > 3 and i < 7):
            iX = 4 + (i - 2) % 3
            iY = 4 + (i ) % 3
        else:
            iX = i
            iY = i
        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[iX][which] = value * e1
                    Y[iY][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))

    # Python indexing starts at zero, so we're always shifted down an
    # array from the L&T numbering scheme:
    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[5]+Y[5]+Z[5], axis=1)
                bc_part = numpy.append(bc_array, X[6]+Y[6]+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:	1917
Committed:	Wed Jul 31 12:58:35 2013 UTC (12 years, 1 month ago) by gezelter
File size:	14806 byte(s)
Log Message:	Fixed all of the dipolar arrays to match the L&T indexing scheme
#	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	gezelter	1917	# The X and Y arrays are cubic rotations of the Z array, so we
77			# need to re-order them to match the numbering scheme in
78			# Luttinger & Tisza:
79			if (i > 0 and i < 4):
80			iX = 1 + (i + 1) % 3
81			iY = 1 + (i ) % 3
82			elif (i > 3 and i < 7):
83			iX = 4 + (i - 2) % 3
84			iY = 4 + (i ) % 3
85			else:
86			iX = i
87			iY = i
88	gezelter	1898	which = 0
89			for l1 in range(2):
90			for l2 in range(2):
91			for l3 in range(2):
92			lvals = numpy.array([l1,l2,l3])
93			value = math.pow(-1, numpy.dot(cell[i], lvals))
94	gezelter	1917	X[iX][which] = value * e1
95			Y[iY][which] = value * e2
96			Z[i][which] = value * e3
97	gezelter	1909	which = which+1
98
99	gezelter	1907
100			lp_array = numpy.zeros(0).reshape((0,3))
101			for i in range(8):
102	gezelter	1909	lp_array = numpy.vstack((lp_array, corners[i]))
103	gezelter	1907
104			bc_array = numpy.zeros(0).reshape((0,3))
105			for i in range(8):
106	gezelter	1909	bc_array = numpy.vstack((bc_array, corners[i] + [0.5,0.5,0.5]))
107	gezelter	1898
108	gezelter	1907	xy_array = numpy.zeros(0).reshape((0,3))
109			for i in range(8):
110	gezelter	1909	xy_array = numpy.vstack((xy_array, corners[i] + [0.5,0.5,0.0]))
111	gezelter	1898
112	gezelter	1907	xz_array = numpy.zeros(0).reshape((0,3))
113			for i in range(8):
114	gezelter	1909	xz_array = numpy.vstack((xz_array, corners[i] + [0.5,0.0,0.5]))
115	gezelter	1898
116	gezelter	1907	yz_array = numpy.zeros(0).reshape((0,3))
117			for i in range(8):
118	gezelter	1909	yz_array = numpy.vstack((yz_array, corners[i] + [0.0,0.5,0.5]))
119	gezelter	1898
120	gezelter	1907	known_case = False
121			basic_array = numpy.zeros(0).reshape((0,3,3))
122	gezelter	1889
123	gezelter	1907	lp_part = numpy.zeros(0).reshape((0,3,3))
124			bc_part = numpy.zeros(0).reshape((0,3,3))
125			xy_part = numpy.zeros(0).reshape((0,3,3))
126			xz_part = numpy.zeros(0).reshape((0,3,3))
127			yz_part = numpy.zeros(0).reshape((0,3,3))
128	gezelter	1889
129	gezelter	1917	# Python indexing starts at zero, so we're always shifted down an
130			# array from the L&T numbering scheme:
131	gezelter	1899	if (arrayType == 'X'):
132			if (int(latticeType)):
133			which = int(latticeType) - 1
134	gezelter	1907	basic_array = numpy.append(lp_array, X[which], axis=1)
135	gezelter	1899	known_case = True
136			if (arrayType == 'Y'):
137			if (int(latticeType)):
138			which = int(latticeType) - 1
139	gezelter	1907	basic_array = numpy.append(lp_array, Y[which], axis=1)
140	gezelter	1899	known_case = True
141			if (arrayType == 'Z'):
142			if (int(latticeType)):
143			which = int(latticeType) - 1
144	gezelter	1907	basic_array = numpy.append(lp_array, Z[which], axis=1)
145	gezelter	1899	known_case = True
146			if (arrayType == 'A'):
147	gezelter	1898	if (latticeType.lower() == 'sc'):
148	gezelter	1907	basic_array = numpy.append(lp_array, Z[4], axis=1)
149	gezelter	1898	known_case = True
150			if (latticeType.lower() == 'bcc'):
151			if (dipoleDirection.lower() == '001'):
152	gezelter	1907	lp_part = numpy.append(lp_array, Z[0], axis=1)
153			bc_part = numpy.append(bc_array, -Z[0], axis=1)
154			basic_array = numpy.append(lp_part, bc_part, axis=0)
155	gezelter	1898	known_case = True
156			if (dipoleDirection.lower() == '111'):
157	gezelter	1907	lp_part = numpy.append(lp_array, Z[4]+X[6]+Y[5], axis=1)
158			bc_part = numpy.append(bc_array, Z[4]+X[6]+Y[5], axis=1)
159			basic_array = numpy.append(lp_part, bc_part, axis=0)
160	gezelter	1898	known_case = True
161	gezelter	1912	if (dipoleDirection.lower() == 'min'):
162	gezelter	1917	lp_part = numpy.append(lp_array, X[5]+Y[5]+Z[5], axis=1)
163			bc_part = numpy.append(bc_array, X[6]+Y[6]+Z[6], axis=1)
164	gezelter	1912	basic_array = numpy.append(lp_part, bc_part, axis=0)
165			known_case = True
166	gezelter	1898	if (latticeType.lower() == 'fcc'):
167			if (dipoleDirection.lower() == '001'):
168	gezelter	1907	lp_part = numpy.append(lp_array, Z[0], axis=1)
169			xy_part = numpy.append(xy_array, Z[0], axis=1)
170			xz_part = numpy.append(xz_array, -Z[0], axis=1)
171			yz_part = numpy.append(yz_array, -Z[0], axis=1)
172
173			basic_array = numpy.append(lp_part, xy_part, axis=0)
174			basic_array = numpy.append(basic_array, xz_part, axis=0)
175			basic_array = numpy.append(basic_array, yz_part, axis=0)
176
177	gezelter	1898	known_case = True
178			if (dipoleDirection.lower() == '011'):
179	gezelter	1907	lp_part = numpy.append(lp_array, Z[0]+Y[0], axis=1)
180			xy_part = numpy.append(xy_array, -Z[0]-Y[0], axis=1)
181			xz_part = numpy.append(xz_array, -Z[0]-Y[0], axis=1)
182			yz_part = numpy.append(yz_array, Z[0]+Y[0], axis=1)
183
184			basic_array = numpy.append(lp_part, xy_part, axis=0)
185			basic_array = numpy.append(basic_array, xz_part, axis=0)
186			basic_array = numpy.append(basic_array, yz_part, axis=0)
187
188	gezelter	1898	known_case = True
189	gezelter	1889	else:
190			if (latticeType.lower() == 'sc'):
191	gezelter	1907	basic_array = numpy.append(lp_array, Z[4], axis=1)
192	gezelter	1898	known_case = True
193			if (latticeType.lower() == 'bcc'):
194	gezelter	1889	if (dipoleDirection.lower() == '001'):
195	gezelter	1907	lp_part = numpy.append(lp_array, Z[4], axis=1)
196			bc_part = numpy.append(bc_array, -Z[4], axis=1)
197			basic_array = numpy.append(lp_part, bc_part, axis=0)
198	gezelter	1898	known_case = True
199			if (dipoleDirection.lower() == '111'):
200	gezelter	1907	lp_part = numpy.append(lp_array, Z[4]+X[6]+Y[5], axis=1)
201			bc_part = numpy.append(bc_array, Z[4]+X[6]+Y[5], axis=1)
202			basic_array = numpy.append(lp_part, bc_part, axis=0)
203	gezelter	1898	known_case = True
204			if (latticeType.lower() == 'fcc'):
205	gezelter	1889	if (dipoleDirection.lower() == '001'):
206	gezelter	1907	lp_part = numpy.append(lp_array, Z[0], axis=1)
207			xy_part = numpy.append(xy_array, -Z[0], axis=1)
208			xz_part = numpy.append(xz_array, -Z[0], axis=1)
209			yz_part = numpy.append(yz_array, Z[0], axis=1)
210
211			basic_array = numpy.append(lp_part, xy_part, axis=0)
212			basic_array = numpy.append(basic_array, xz_part, axis=0)
213			basic_array = numpy.append(basic_array, yz_part, axis=0)
214
215	gezelter	1898	known_case = True
216			if (dipoleDirection.lower() == '011'):
217	gezelter	1907	lp_part = numpy.append(lp_array, Z[7]+Y[7], axis=1)
218			xy_part = numpy.append(xy_array, Z[7]+Y[7], axis=1)
219			xz_part = numpy.append(xz_array, -Z[7]-Y[7], axis=1)
220			yz_part = numpy.append(yz_array, Z[7]+Y[7], axis=1)
221
222			basic_array = numpy.append(lp_part, xy_part, axis=0)
223			basic_array = numpy.append(basic_array, xz_part, axis=0)
224			basic_array = numpy.append(basic_array, yz_part, axis=0)
225
226	gezelter	1898	known_case = True
227	gezelter	1889
228	gezelter	1898	if (not known_case):
229			print "unhandled combination of lattice and dipole direction"
230			print __doc__
231
232	gezelter	1907	bravais_lattice = numpy.zeros(0).reshape((0,6))
233	gezelter	1889	for i in range(latticeNumber):
234			for j in range(latticeNumber):
235			for k in range(latticeNumber):
236	gezelter	1898	for l in range(len(basic_array)):
237	gezelter	1907	lat_vec = numpy.array([[2i, 2j, 2*k, 0.0, 0.0, 0.0]])
238			bravais_lattice = numpy.append(bravais_lattice, lat_vec + basic_array[l], axis=0)
239	gezelter	1889
240	gezelter	1898	outputFile = open(outputFileName, 'w')
241
242			outputFile.write('<OpenMD version=2>\n')
243			outputFile.write(' <MetaData>\n')
244			outputFile.write(' molecule{\n')
245			outputFile.write(' name = \"D\";\n')
246			outputFile.write(' \n')
247			outputFile.write(' atom[0]{\n')
248			outputFile.write(' type = \"D\";\n')
249			outputFile.write(' position(0.0, 0.0, 0.0);\n')
250			outputFile.write(' orientation(0.0, 0.0, 0.0);\n')
251			outputFile.write(' }\n')
252			outputFile.write(' }\n')
253			outputFile.write(' component{\n')
254			outputFile.write(' type = \"D\";\n')
255			outputFile.write(' nMol = '+ repr(len(bravais_lattice)) + ';\n')
256			outputFile.write(' }\n')
257	gezelter	1889
258	gezelter	1898	outputFile.write(' ensemble = NVE;\n')
259			outputFile.write(' forceField = \"Multipole\";\n')
260
261			outputFile.write(' cutoffMethod = \"shifted_force\";\n')
262			outputFile.write(' electrostaticScreeningMethod = \"damped\";\n')
263
264			outputFile.write(' cutoffRadius = 9.0;\n')
265			outputFile.write(' dampingAlpha = 0.18;\n')
266			outputFile.write(' statFileFormat = \"TIME\|TOTAL_ENERGY\|POTENTIAL_ENERGY\|KINETIC_ENERGY\|TEMPERATURE\|PRESSURE\|VOLUME\|CONSERVED_QUANTITY\|ELECTROSTATIC_POTENTIAL\";\n')
267			outputFile.write(' dt = 1.0;\n')
268			outputFile.write(' runTime = 1.0;\n')
269			outputFile.write(' sampleTime = 1.0;\n')
270			outputFile.write(' statusTime = 1.0;\n')
271			outputFile.write(' </MetaData>\n')
272			outputFile.write(' <Snapshot>\n')
273			outputFile.write(' <FrameData>\n');
274			outputFile.write(" Time: %.10g\n" % (0.0))
275	gezelter	1889
276	gezelter	1898	Hxx = 2.0 * latticeConstant * latticeNumber
277			Hyy = 2.0 * latticeConstant * latticeNumber
278			Hzz = 2.0 * latticeConstant * latticeNumber
279
280			outputFile.write(' Hmat: {{%d, 0, 0}, {0, %d, 0}, {0, 0, %d}}\n' % (Hxx, Hyy, Hzz))
281			outputFile.write(' </FrameData>\n')
282			outputFile.write(' <StuntDoubles>\n')
283			sdFormat = 'pvqj'
284			index = 0
285	gezelter	1907
286	gezelter	1889	for i in range(len(bravais_lattice)):
287			xcart = latticeConstant*(bravais_lattice[i][0])
288			ycart = latticeConstant*(bravais_lattice[i][1])
289			zcart = latticeConstant*(bravais_lattice[i][2])
290	gezelter	1898	dx = bravais_lattice[i][3]
291			dy = bravais_lattice[i][4]
292			dz = bravais_lattice[i][5]
293
294	gezelter	1909	dlen = math.sqrt(dxdx + dydy + dz*dz)
295			ctheta = dz / dlen
296			theta = math.acos(ctheta)
297			stheta = math.sqrt(1.0 - ctheta*ctheta)
298			psi = 0.0
299	gezelter	1914	phi = math.atan2(dx/dlen, -dy/dlen)
300	gezelter	1898
301	gezelter	1909	q = [0.0,0.0,0.0,0.0]
302			q[0] = math.cos(theta/2)*math.cos((phi+psi)/2)
303			q[1] = math.sin(theta/2)*math.cos((phi-psi)/2)
304			q[2] = math.sin(theta/2)*math.sin((phi-psi)/2)
305			q[3] = math.cos(theta/2)*math.sin((phi+psi)/2)
306	gezelter	1898
307	gezelter	1905	qlen = math.sqrt(q[0]q[0] + q[1]q[1] + q[2]q[2] + q[3]q[3])
308			q[0] = q[0]/qlen
309			q[1] = q[1]/qlen
310			q[2] = q[2]/qlen
311			q[3] = q[3]/qlen
312	gezelter	1898
313	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))
314	gezelter	1898	index = index+1
315
316			outputFile.write(" </StuntDoubles>\n")
317			outputFile.write(" </Snapshot>\n")
318			outputFile.write("</OpenMD>\n")
319	gezelter	1889	outputFile.close()
320
321	gezelter	1898	outputFile.close()
322
323	gezelter	1889	def main(argv):
324
325	gezelter	1899	arrayType = "A"
326	gezelter	1889	haveOutputFileName = False
327			latticeType = "fcc"
328			dipoleDirection = "001"
329			latticeNumber = 3
330			latticeConstant = 4
331			try:
332	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="])
333	gezelter	1889	except getopt.GetoptError:
334			usage()
335			sys.exit(2)
336			for opt, arg in opts:
337			if opt in ("-h", "--help"):
338			usage()
339			sys.exit()
340	gezelter	1899	elif opt in ("-x", "--array-type-X"):
341			arrayType = "X"
342			elif opt in ("-y", "--array-type-Y"):
343			arrayType = "Y"
344			elif opt in ("-z", "--array-type-Z"):
345			arrayType = "Z"
346	gezelter	1889	elif opt in ("-b", "--array-type-B"):
347	gezelter	1899	arrayType = "B"
348	gezelter	1889	elif opt in ("-l", "--lattice"):
349			latticeType = arg
350			elif opt in ("-d", "--direction"):
351			dipoleDirection = arg
352			elif opt in ("-c", "--constant"):
353			latticeConstant = float(arg)
354			elif opt in ("-n"):
355			latticeNumber = int(arg)
356			elif opt in ("-o", "--output-file"):
357			outputFileName = arg
358			haveOutputFileName = True
359			if (not haveOutputFileName):
360			usage()
361			print "No output file was specified"
362			sys.exit()
363
364	gezelter	1899	createLattice(latticeType, latticeNumber, latticeConstant, dipoleDirection, arrayType, outputFileName);
365	gezelter	1889
366			if __name__ == "__main__":
367			if len(sys.argv) == 1:
368			usage()
369			sys.exit()
370			main(sys.argv[1:])