Révision 12 ase/calculators/mopac.py
| mopac.py (revision 12) | ||
|---|---|---|
| 8 | 8 |
import commands |
| 9 | 9 |
import numpy as np |
| 10 | 10 |
|
| 11 |
from copy import deepcopy |
|
| 12 |
|
|
| 11 | 13 |
from ase.atoms import Atoms |
| 12 | 14 |
from ase.units import kcal, mol |
| 13 | 15 |
from ase.calculators.general import Calculator |
| ... | ... | |
| 22 | 24 |
self.spin = 0 |
| 23 | 25 |
self.occupations = None |
| 24 | 26 |
|
| 25 |
self.first = True |
|
| 26 | 27 |
self.stress = None |
| 27 | 28 |
self.energy_zero = None |
| 28 | 29 |
self.energy_free = None |
| 29 | 30 |
self.forces = None |
| 30 | 31 |
|
| 32 |
self.atoms_old = None |
|
| 33 |
|
|
| 31 | 34 |
def write_moldata(self): |
| 32 | 35 |
""" |
| 33 | 36 |
Writes the files that have to be written each timestep |
| ... | ... | |
| 47 | 50 |
|
| 48 | 51 |
def run(self): |
| 49 | 52 |
self.write_moldata() |
| 53 |
print '<MOPAC calculation started>' |
|
| 50 | 54 |
out = commands.getoutput('mopac pymopac.in')
|
| 51 |
|
|
| 52 | 55 |
outputfile = open('pymopac.in.out')
|
| 53 | 56 |
outfile = open('mopac%02i.out'%(self.restart),'w')
|
| 54 | 57 |
for line in outputfile: |
| ... | ... | |
| 60 | 63 |
if energy == None: |
| 61 | 64 |
energy = self.rerun() |
| 62 | 65 |
self.forces = self.read_forces() |
| 66 |
print '<MOPAC calculation finished>' |
|
| 63 | 67 |
|
| 64 | 68 |
self.energy_zero = energy |
| 65 | 69 |
self.energy_free = energy |
| ... | ... | |
| 98 | 102 |
return energy |
| 99 | 103 |
|
| 100 | 104 |
def read_forces(self): |
| 101 |
|
|
| 102 | 105 |
outfile = open('pymopac.in.aux')
|
| 103 | 106 |
lines = outfile.readlines() |
| 104 | 107 |
outfile.close() |
| 105 | 108 |
|
| 106 | 109 |
outputforces = None |
| 107 | 110 |
nats = len(self.atoms) |
| 108 |
nx = nats * 3 |
|
| 109 |
|
|
| 111 |
|
|
| 112 |
forces = np.zeros((nats, 3), float) |
|
| 110 | 113 |
for i,line in enumerate(lines): |
| 111 | 114 |
if line.find('GRADIENTS:') != -1:
|
| 112 |
forces = [] |
|
| 113 | 115 |
l = 0 |
| 114 |
for j in range(nx):
|
|
| 116 |
for j in range(nats * 3):
|
|
| 115 | 117 |
if j%10 == 0: |
| 116 | 118 |
l += 1 |
| 117 | 119 |
gline = lines[i+l] |
| 118 | 120 |
k = j -((l-1)*10) |
| 121 |
|
|
| 119 | 122 |
try: |
| 120 |
forces.append(-float(gline[k*18:(k*18)+18]))
|
|
| 123 |
f = -float(gline[k*18:(k*18)+18])
|
|
| 121 | 124 |
except ValueError: |
| 122 | 125 |
if outputforces == None: |
| 123 | 126 |
outputforces = self.read_output_forces() |
| 124 |
forces.append(outputforces[j]) |
|
| 127 |
f = outputforces[j] |
|
| 128 |
forces[j/3, j%3] = f |
|
| 125 | 129 |
break |
| 126 | 130 |
|
| 127 |
forces = np.reshape(forces, (nats, 3)) |
|
| 128 | 131 |
forces *= kcal_mol |
| 129 | 132 |
return forces |
| 130 | 133 |
|
| ... | ... | |
| 134 | 137 |
outfile.close() |
| 135 | 138 |
|
| 136 | 139 |
nats = len(self.atoms) |
| 137 |
nx = nats * 3 |
|
| 138 | 140 |
|
| 139 |
forces = []
|
|
| 141 |
forces = np.zeros((nats, 3), float)
|
|
| 140 | 142 |
for i,line in enumerate(lines): |
| 141 | 143 |
if line.find('GRADIENT\n') != -1:
|
| 142 |
for j in range(nx):
|
|
| 144 |
for j in range(nats * 3):
|
|
| 143 | 145 |
gline = lines[i+j+1] |
| 144 |
forces.append(-float(gline[49:62]))
|
|
| 146 |
forces[j/3, j%3] = -float(gline[49:62])
|
|
| 145 | 147 |
break |
| 146 | 148 |
|
| 147 | 149 |
return forces |
| ... | ... | |
| 183 | 185 |
return block |
| 184 | 186 |
|
| 185 | 187 |
def update(self, atoms): |
| 186 |
if not (self.atoms.get_positions() == atoms.get_positions()).all() or self.first:
|
|
| 188 |
if self.atoms_old != atoms:
|
|
| 187 | 189 |
self.atoms = atoms.copy() |
| 190 |
self.atoms_old = atoms.copy() |
|
| 188 | 191 |
self.run() |
Formats disponibles : Unified diff