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