modified code to accommodate partial charges; devel01 in #15

PyCT/core.py

@@ -95,7 +95,7 @@ def __init__(self, materialParameters):
         self.name = materialParameters.name
         self.speciesTypes = materialParameters.speciesTypes[:]
         self.numSpeciesTypes = len(self.speciesTypes)
-        self.speciesChargeList = materialParameters.speciesChargeList[:]
+        self.speciesChargeList = deepcopy(materialParameters.speciesChargeList)
         self.speciesToElementTypeMap = deepcopy(
                                     materialParameters.speciesToElementTypeMap)
 
@@ -953,21 +953,11 @@ def __init__(self, materialInfo, materialNeighbors, hopNeighborList,
 
         self.pbc = self.neighbors.pbc
         self.speciesCount = speciesCount
-        self.systemCharge = np.dot(speciesCount,
-                                   self.material.speciesChargeList)
 
         # total number of unit cells
         self.systemSize = self.neighbors.systemSize
         self.numCells = self.systemSize.prod()
 
-        # generate lattice charge list
-        unitcellChargeList = np.array(
-                [self.material.chargeTypes[
-                    self.material.elementTypes[elementTypeIndex]]
-                 for elementTypeIndex in self.material.elementTypeIndexList])
-        self.latticeChargeList = np.tile(
-                            unitcellChargeList, self.numCells)[:, np.newaxis]
-
         self.cumulativeDisplacementList = cumulativeDisplacementList
 
         # variables for ewald sum
@@ -1028,17 +1018,24 @@ def generateRandomOccupancy(self, speciesCount):
                                         numSpecies)[:])
         return occupancy
 
-    def chargeConfig(self, occupancy):
+    def chargeConfig(self, occupancy, ionChargeType, speciesChargeType):
         """Returns charge distribution of the current configuration"""
-        chargeList = np.copy(self.latticeChargeList)
+
+        # generate lattice charge list
+        unitcellChargeList = np.array(
+                [self.material.chargeTypes[ionChargeType][
+                    self.material.elementTypes[elementTypeIndex]]
+                 for elementTypeIndex in self.material.elementTypeIndexList])
+        chargeList = np.tile(unitcellChargeList, self.numCells)[:, np.newaxis]
 
         for speciesTypeIndex in range(self.material.numSpeciesTypes):
             startIndex = 0 + self.speciesCount[:speciesTypeIndex].sum()
             endIndex = startIndex + self.speciesCount[speciesTypeIndex]
             centerSiteSystemElementIndices = occupancy[
                                                     startIndex:endIndex][:]
-            chargeList[centerSiteSystemElementIndices] += \
-                self.material.speciesChargeList[speciesTypeIndex]
+            chargeList[centerSiteSystemElementIndices] += (
+                        self.material.speciesChargeList[speciesChargeType][
+                                                            speciesTypeIndex])
         return chargeList
 
     def ewaldSumSetup(self, outdir=None):
@@ -1110,8 +1107,8 @@ def generatePreComputedArrayLogReport(self, outdir):
 class run(object):
     """defines the subroutines for running Kinetic Monte Carlo and
         computing electrostatic interaction energies"""
-    def __init__(self, system, precomputedArray, T, nTraj,
-                 tFinal, timeInterval):
+    def __init__(self, system, precomputedArray, T, ionChargeType,
+                 speciesChargeType, nTraj, tFinal, timeInterval):
         """Returns the PBC condition of the system"""
         self.startTime = datetime.now()
 
@@ -1120,6 +1117,8 @@ def __init__(self, system, precomputedArray, T, nTraj,
         self.neighbors = self.system.neighbors
         self.precomputedArray = precomputedArray
         self.T = T * self.material.K2AUTEMP
+        self.ionChargeType = ionChargeType
+        self.speciesChargeType = speciesChargeType
         self.nTraj = int(nTraj)
         self.tFinal = tFinal * self.material.SEC2AUTIME
         self.timeInterval = timeInterval * self.material.SEC2AUTIME
@@ -1140,8 +1139,9 @@ def __init__(self, system, precomputedArray, T, nTraj,
                         index
                         for index, value in enumerate(self.system.speciesCount)
                         for _ in range(value)]
-        self.speciesChargeList = [self.material.speciesChargeList[index]
-                                  for index in self.speciesTypeIndexList]
+        self.speciesChargeList = [
+                self.material.speciesChargeList[self.speciesChargeType][index]
+                for index in self.speciesTypeIndexList]
         self.hopElementTypeList = [
                                 self.material.hopElementTypes[speciesType][0]
                                 for speciesType in self.speciesTypeList]
@@ -1235,16 +1235,21 @@ def doKMCSteps(self, outdir, report=1, randomSeed=1):
         neighborSiteSystemElementIndexList = np.zeros(self.nProc, dtype=int)
         rowIndexList = np.zeros(self.nProc, dtype=int)
         neighborIndexList = np.zeros(self.nProc, dtype=int)
+        systemCharge = np.dot(
+                    self.system.speciesCount,
+                    self.material.speciesChargeList[self.speciesChargeType])
 
         ewaldNeut = - (np.pi
-                       * (self.system.systemCharge**2)
+                       * (systemCharge**2)
                        / (2 * self.system.systemVolume * self.system.alpha))
         precomputedArray = self.precomputedArray
         for _ in range(nTraj):
             currentStateOccupancy = self.system.generateRandomOccupancy(
                                                     self.system.speciesCount)
             currentStateChargeConfig = self.system.chargeConfig(
-                                                        currentStateOccupancy)
+                                                        currentStateOccupancy,
+                                                        self.ionChargeType,
+                                                        self.speciesChargeType)
             currentStateChargeConfigProd = np.multiply(
                                         currentStateChargeConfig.transpose(),
                                         currentStateChargeConfig)

PyCT/materialMSD.py

@@ -8,8 +8,9 @@
 
 
 def materialMSD(systemDirectoryPath, fileFormatIndex, systemSize, pbc, nDim,
-                Temp, speciesCount, tFinal, nTraj, timeInterval, msdTFinal,
-                trimLength, displayErrorBars, reprTime, reprDist, report):
+                Temp, ionChargeType, speciesChargeType, speciesCount, tFinal,
+                nTraj, timeInterval, msdTFinal, trimLength, displayErrorBars,
+                reprTime, reprDist, report):
 
     # Load material parameters
     configDirName = 'ConfigurationFiles'
@@ -34,16 +35,18 @@ def materialMSD(systemDirectoryPath, fileFormatIndex, systemSize, pbc, nDim,
 
     # Change to working directory
     parentDir1 = 'SimulationFiles'
+    parentDir2 = ('ionChargeType=' + ionChargeType
+                  + '; speciesChargeType=' + speciesChargeType)
     nElectrons = speciesCount[0]
     nHoles = speciesCount[1]
-    parentDir2 = (str(nElectrons)
+    parentDir3 = (str(nElectrons)
                   + ('electron' if nElectrons == 1 else 'electrons') + ', '
                   + str(nHoles) + ('hole' if nHoles == 1 else 'holes'))
-    parentDir3 = str(Temp) + 'K'
+    parentDir4 = str(Temp) + 'K'
     workDir = (('%1.2E' % tFinal) + 'SEC,' + ('%1.2E' % timeInterval)
                + 'TimeInterval,' + ('%1.2E' % nTraj) + 'Traj')
     workDirPath = os.path.join(systemDirectoryPath, parentDir1, parentDir2,
-                               parentDir3, workDir)
+                               parentDir3, parentDir4, workDir)
 
     if not os.path.exists(workDirPath):
         print('Simulation files do not exist. Aborting.')

PyCT/materialRun.py

@@ -9,8 +9,8 @@
 
 
 def materialRun(systemDirectoryPath, fileFormatIndex, systemSize, pbc, Temp,
-                speciesCount, tFinal, nTraj, timeInterval, randomSeed, report,
-                overWrite):
+                ionChargeType, speciesChargeType, speciesCount, tFinal, nTraj,
+                timeInterval, randomSeed, report, overWrite):
 
     # Load material parameters
     configDirName = 'ConfigurationFiles'
@@ -38,16 +38,18 @@ def materialRun(systemDirectoryPath, fileFormatIndex, systemSize, pbc, Temp,
 
     # Change to working directory
     parentDir1 = 'SimulationFiles'
+    parentDir2 = ('ionChargeType=' + ionChargeType
+                  + '; speciesChargeType=' + speciesChargeType)
     nElectrons = speciesCount[0]
     nHoles = speciesCount[1]
-    parentDir2 = (str(nElectrons)
+    parentDir3 = (str(nElectrons)
                   + ('electron' if nElectrons == 1 else 'electrons') + ', '
                   + str(nHoles) + ('hole' if nHoles == 1 else 'holes'))
-    parentDir3 = str(Temp) + 'K'
+    parentDir4 = str(Temp) + 'K'
     workDir = (('%1.2E' % tFinal) + 'SEC,' + ('%1.2E' % timeInterval)
                + 'TimeInterval,' + ('%1.2E' % nTraj) + 'Traj')
     workDirPath = os.path.join(systemDirectoryPath, parentDir1, parentDir2,
-                               parentDir3, workDir)
+                               parentDir3, parentDir4, workDir)
     if not os.path.exists(workDirPath):
         os.makedirs(workDirPath)
     os.chdir(workDirPath)
@@ -85,8 +87,9 @@ def materialRun(systemDirectoryPath, fileFormatIndex, systemSize, pbc, Temp,
         precomputedArrayFilePath = os.path.join(inputFileDirectoryPath,
                                                 'precomputedArray.npy')
         precomputedArray = np.load(precomputedArrayFilePath)
-        materialRun = run(materialSystem, precomputedArray, Temp, nTraj,
-                          tFinal, timeInterval)
+        materialRun = run(materialSystem, precomputedArray, Temp,
+                          ionChargeType, speciesChargeType, nTraj, tFinal,
+                          timeInterval)
 
         materialRun.doKMCSteps(workDirPath, report, randomSeed)
     else: