[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__ = "$Revision: 1639 $"
__date__ = "$Date: 2011-09-24 16:18:07 -0400 (Sat, 24 Sep 2011) $"

__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]])

    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, cell[i]))
    
    bc_array = numpy.zeros(0).reshape((0,3))
    for i in range(8):
        bc_array = numpy.vstack((bc_array, cell[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, cell[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, cell[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, cell[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

    print bravais_lattice
    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]

        uz = numpy.array([dx, dy, dz])
        uz = uz/numpy.linalg.norm(uz)
        
        uy = numpy.array([0.0, 1.0, 0.0])
        uy = uy - uz * numpy.vdot(uy, uz) / numpy.vdot(uz, uz) 
        uy = uy/numpy.linalg.norm(uy)

        ux = numpy.cross(uy, uz)

        RotMat = [ux, uy, uz]

        q = [0.0, 0.0, 0.0, 0.0]

        # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:

        t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
        
        if( t > 1e-6 ):
            s = 0.5 / math.sqrt( t )
            q[0] = 0.25 / s
            q[1] = (RotMat[1][2] - RotMat[2][1]) * s
            q[2] = (RotMat[2][0] - RotMat[0][2]) * s
            q[3] = (RotMat[0][1] - RotMat[1][0]) * s
        else:
            ad1 = RotMat[0][0]
            ad2 = RotMat[1][1]
            ad3 = RotMat[2][2]

            if( ad1 >= ad2 and ad1 >= ad3 ):
                s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
                q[0] = (RotMat[1][2] - RotMat[2][1]) * s
                q[1] = 0.25 / s
                q[2] = (RotMat[0][1] + RotMat[1][0]) * s
                q[3] = (RotMat[0][2] + RotMat[2][0]) * s
            elif ( ad2 >= ad1 and ad2 >= ad3 ):
                s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
                q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
                q[1] = (RotMat[0][1] + RotMat[1][0]) * s
                q[2] = 0.25 / s
                q[3] = (RotMat[1][2] + RotMat[2][1]) * s
            else:
                s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
                q[0] = (RotMat[0][1] - RotMat[1][0]) * s
                q[1] = (RotMat[0][2] + RotMat[2][0]) * s
                q[2] = (RotMat[1][2] + RotMat[2][1]) * s
                q[3] = 0.25 / s

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