Révision 3 ase/calculators/turbomole.py
| turbomole.py (revision 3) | ||
|---|---|---|
| 2 | 2 |
|
| 3 | 3 |
http://www.turbomole.com/ |
| 4 | 4 |
""" |
| 5 |
|
|
| 6 |
|
|
| 7 | 5 |
import os, sys, string |
| 8 | 6 |
|
| 9 |
from ase.data import chemical_symbols |
|
| 10 | 7 |
from ase.units import Hartree, Bohr |
| 11 | 8 |
from ase.io.turbomole import write_turbomole |
| 12 |
from general import Calculator |
|
| 9 |
from ase.calculators.general import Calculator
|
|
| 13 | 10 |
|
| 14 | 11 |
import numpy as np |
| 15 | 12 |
|
| 16 | 13 |
class Turbomole(Calculator): |
| 17 | 14 |
"""Class for doing Turbomole calculations. |
| 18 | 15 |
""" |
| 19 |
def __init__(self, label='turbomole'): |
|
| 16 |
def __init__(self, label='turbomole', calculate_energy="dscf", calculate_forces="grad", post_HF = False):
|
|
| 20 | 17 |
"""Construct TURBOMOLE-calculator object. |
| 21 | 18 |
|
| 22 | 19 |
Parameters |
| ... | ... | |
| 59 | 56 |
e = a.get_potential_energy() |
| 60 | 57 |
|
| 61 | 58 |
""" |
| 62 |
|
|
| 63 |
# get names of executables for turbomole energy and forces |
|
| 64 |
# get the path |
|
| 65 |
|
|
| 66 |
os.system('rm -f sysname.file; sysname > sysname.file')
|
|
| 67 |
f = open('sysname.file')
|
|
| 68 |
architechture = f.readline()[:-1] |
|
| 69 |
f.close() |
|
| 70 |
tmpath = os.environ['TURBODIR'] |
|
| 71 |
pre = tmpath + '/bin/' + architechture + '/' |
|
| 72 |
|
|
| 73 |
|
|
| 74 |
if os.path.isfile('control'):
|
|
| 75 |
f = open('control')
|
|
| 76 |
else: |
|
| 77 |
print 'File control is missing' |
|
| 78 |
raise RuntimeError, \ |
|
| 79 |
'Please run Turbomole define and come thereafter back' |
|
| 80 |
lines = f.readlines() |
|
| 81 |
f.close() |
|
| 82 |
self.tm_program_energy=pre+'dscf' |
|
| 83 |
self.tm_program_forces=pre+'grad' |
|
| 84 |
for line in lines: |
|
| 85 |
if line.startswith('$ricore'):
|
|
| 86 |
self.tm_program_energy=pre+'ridft' |
|
| 87 |
self.tm_program_forces=pre+'rdgrad' |
|
| 88 |
|
|
| 59 |
|
|
| 89 | 60 |
self.label = label |
| 90 | 61 |
self.converged = False |
| 91 |
#clean up turbomole energy file |
|
| 92 |
os.system('rm -f energy; touch energy')
|
|
| 62 |
|
|
| 63 |
# set calculators for energy and forces |
|
| 64 |
self.calculate_energy = calculate_energy |
|
| 65 |
self.calculate_forces = calculate_forces |
|
| 93 | 66 |
|
| 94 |
#turbomole has no stress |
|
| 67 |
# turbomole has no stress
|
|
| 95 | 68 |
self.stress = np.empty((3, 3)) |
| 96 | 69 |
|
| 70 |
# storage for energy and forces |
|
| 71 |
self.e_total = None |
|
| 72 |
self.forces = None |
|
| 73 |
self.updated = False |
|
| 74 |
|
|
| 75 |
# atoms must be set |
|
| 76 |
self.atoms = None |
|
| 77 |
|
|
| 78 |
# POST-HF method |
|
| 79 |
self.post_HF = post_HF |
|
| 97 | 80 |
|
| 98 |
def update(self, atoms): |
|
| 99 |
"""Energy and forces are calculated when atoms have moved |
|
| 100 |
by calling self.calculate |
|
| 101 |
""" |
|
| 102 |
if (not self.converged or |
|
| 103 |
len(self.numbers) != len(atoms) or |
|
| 104 |
(self.numbers != atoms.get_atomic_numbers()).any()): |
|
| 105 |
self.initialize(atoms) |
|
| 106 |
self.calculate(atoms) |
|
| 107 |
elif ((self.positions != atoms.get_positions()).any() or |
|
| 108 |
(self.pbc != atoms.get_pbc()).any() or |
|
| 109 |
(self.cell != atoms.get_cell()).any()): |
|
| 110 |
self.calculate(atoms) |
|
| 111 |
|
|
| 112 | 81 |
def initialize(self, atoms): |
| 113 | 82 |
self.numbers = atoms.get_atomic_numbers().copy() |
| 114 | 83 |
self.species = [] |
| ... | ... | |
| 116 | 85 |
self.species.append(Z) |
| 117 | 86 |
self.converged = False |
| 118 | 87 |
|
| 88 |
def execute(self, command): |
|
| 89 |
from subprocess import Popen, PIPE |
|
| 90 |
try: |
|
| 91 |
proc = Popen([command], shell=True, stderr=PIPE) |
|
| 92 |
error = proc.communicate()[1] |
|
| 93 |
if "abnormally" in error: |
|
| 94 |
raise OSError(error) |
|
| 95 |
print "TM command: ", command, "successfully executed" |
|
| 96 |
except OSError, e: |
|
| 97 |
print >>sys.stderr, "Execution failed:", e |
|
| 98 |
sys.exit(1) |
|
| 99 |
|
|
| 119 | 100 |
def get_potential_energy(self, atoms): |
| 120 |
self.update(atoms) |
|
| 121 |
return self.etotal |
|
| 101 |
# update atoms |
|
| 102 |
self.set_atoms(atoms) |
|
| 103 |
# if update of energy is neccessary |
|
| 104 |
if self.update_energy: |
|
| 105 |
# calculate energy |
|
| 106 |
self.execute(self.calculate_energy + " > ASE.TM.energy.out") |
|
| 107 |
# check for convergence of dscf cycle |
|
| 108 |
if os.path.isfile('dscf_problem'):
|
|
| 109 |
print 'Turbomole scf energy calculation did not converge' |
|
| 110 |
raise RuntimeError, \ |
|
| 111 |
'Please run Turbomole define and come thereafter back' |
|
| 112 |
# read energy |
|
| 113 |
self.read_energy() |
|
| 114 |
else : |
|
| 115 |
print "taking old values (E)" |
|
| 116 |
self.update_energy = False |
|
| 117 |
return self.e_total |
|
| 122 | 118 |
|
| 123 | 119 |
def get_forces(self, atoms): |
| 124 |
self.update(atoms) |
|
| 120 |
# update atoms |
|
| 121 |
self.set_atoms(atoms) |
|
| 122 |
# complete energy calculations |
|
| 123 |
if self.update_energy: |
|
| 124 |
self.get_potential_energy(atoms) |
|
| 125 |
# if update of forces is neccessary |
|
| 126 |
if self.update_forces: |
|
| 127 |
# calculate forces |
|
| 128 |
self.execute(self.calculate_forces + " > ASE.TM.forces.out") |
|
| 129 |
# read forces |
|
| 130 |
self.read_forces() |
|
| 131 |
else : |
|
| 132 |
print "taking old values (F)" |
|
| 133 |
self.update_forces = False |
|
| 125 | 134 |
return self.forces.copy() |
| 126 | 135 |
|
| 127 | 136 |
def get_stress(self, atoms): |
| 128 |
self.update(atoms) |
|
| 129 | 137 |
return self.stress.copy() |
| 130 |
|
|
| 131 |
def calculate(self, atoms): |
|
| 132 |
"""Total Turbomole energy is calculated (to file 'energy' |
|
| 133 |
also forces are calculated (to file 'gradient') |
|
| 134 |
""" |
|
| 135 |
self.positions = atoms.get_positions().copy() |
|
| 136 |
self.cell = atoms.get_cell().copy() |
|
| 137 |
self.pbc = atoms.get_pbc().copy() |
|
| 138 |
|
|
| 139 |
#write current coordinates to file 'coord' for Turbomole |
|
| 138 |
|
|
| 139 |
def set_atoms(self, atoms): |
|
| 140 |
if self.atoms == atoms: |
|
| 141 |
return |
|
| 142 |
# performs an update of the atoms |
|
| 143 |
super(Turbomole, self).set_atoms(atoms) |
|
| 140 | 144 |
write_turbomole('coord', atoms)
|
| 141 |
|
|
| 142 |
|
|
| 143 |
#Turbomole energy calculation |
|
| 144 |
os.system('rm -f output.energy.dummy; \
|
|
| 145 |
'+ self.tm_program_energy +'> output.energy.dummy') |
|
| 146 |
|
|
| 147 |
#check that the energy run converged |
|
| 148 |
if os.path.isfile('dscf_problem'):
|
|
| 149 |
print 'Turbomole scf energy calculation did not converge' |
|
| 150 |
print 'issue command t2x -c > last.xyz' |
|
| 151 |
print 'and check geometry last.xyz and job.xxx or statistics' |
|
| 152 |
raise RuntimeError, \ |
|
| 153 |
'Please run Turbomole define and come thereafter back' |
|
| 154 |
|
|
| 155 |
|
|
| 156 |
self.read_energy() |
|
| 157 |
|
|
| 158 |
#Turbomole atomic forces calculation |
|
| 159 |
#killing the previous gradient file because |
|
| 160 |
#turbomole gradients are affected by the previous values |
|
| 161 |
os.system('rm -f gradient; rm -f output.forces.dummy; \
|
|
| 162 |
'+ self.tm_program_forces +'> output.forces.dummy') |
|
| 163 |
|
|
| 164 |
self.read_forces(atoms) |
|
| 165 |
|
|
| 166 |
self.converged = True |
|
| 167 |
|
|
| 145 |
# energy and forces must be re-calculated |
|
| 146 |
self.update_energy = True |
|
| 147 |
self.update_forces = True |
|
| 168 | 148 |
|
| 169 | 149 |
def read_energy(self): |
| 170 | 150 |
"""Read Energy from Turbomole energy file.""" |
| ... | ... | |
| 178 | 158 |
elif line.startswith('$'):
|
| 179 | 159 |
pass |
| 180 | 160 |
else: |
| 181 |
#print 'LINE',line |
|
| 161 |
# print 'LINE',line
|
|
| 182 | 162 |
energy_tmp = float(line.split()[1]) |
| 183 |
#print 'energy_tmp',energy_tmp |
|
| 184 |
self.etotal = energy_tmp * Hartree |
|
| 163 |
if self.post_HF: |
|
| 164 |
energy_tmp += float(line.split()[4]) |
|
| 165 |
self.e_total = energy_tmp * Hartree |
|
| 166 |
|
|
| 185 | 167 |
|
| 186 |
|
|
| 187 |
def read_forces(self,atoms): |
|
| 168 |
def read_forces(self): |
|
| 188 | 169 |
"""Read Forces from Turbomole gradient file.""" |
| 189 |
|
|
| 190 | 170 |
file = open('gradient','r')
|
| 191 |
line=file.readline() |
|
| 192 |
line=file.readline() |
|
| 193 |
tmpforces = np.array([[0, 0, 0]]) |
|
| 194 |
while line: |
|
| 195 |
if 'cycle' in line: |
|
| 196 |
for i, dummy in enumerate(atoms): |
|
| 197 |
line=file.readline() |
|
| 198 |
forces_tmp=[] |
|
| 199 |
for i, dummy in enumerate(atoms): |
|
| 200 |
line=file.readline() |
|
| 201 |
line2=string.replace(line,'D','E') |
|
| 202 |
#tmp=np.append(forces_tmp,np.array\ |
|
| 203 |
# ([[float(f) for f in line2.split()[0:3]]])) |
|
| 204 |
tmp=np.array\ |
|
| 205 |
([[float(f) for f in line2.split()[0:3]]]) |
|
| 206 |
tmpforces=np.concatenate((tmpforces,tmp)) |
|
| 207 |
line=file.readline() |
|
| 208 |
|
|
| 171 |
lines=file.readlines() |
|
| 172 |
file.close() |
|
| 209 | 173 |
|
| 210 |
#note the '-' sign for turbomole, to get forces |
|
| 211 |
self.forces = (-np.delete(tmpforces, np.s_[0:1], axis=0))*Hartree/Bohr |
|
| 174 |
forces = np.array([[0, 0, 0]]) |
|
| 175 |
|
|
| 176 |
nline = len(lines) |
|
| 177 |
iline = -1 |
|
| 178 |
|
|
| 179 |
for i in xrange(nline): |
|
| 180 |
if 'cycle' in lines[i]: |
|
| 181 |
iline = i |
|
| 182 |
|
|
| 183 |
if iline < 0: |
|
| 184 |
print 'Gradients not found' |
|
| 185 |
raise RuntimeError("Please check TURBOMOLE gradients")
|
|
| 212 | 186 |
|
| 213 |
#print 'forces', self.forces |
|
| 214 |
|
|
| 215 |
def read(self): |
|
| 216 |
"""Dummy stress for turbomole""" |
|
| 217 |
self.stress = np.empty((3, 3)) |
|
| 187 |
iline += len(self.atoms) + 1 |
|
| 188 |
nline -= 1 |
|
| 189 |
for i in xrange(iline, nline): |
|
| 190 |
line = string.replace(lines[i],'D','E') |
|
| 191 |
#tmp=np.append(forces_tmp,np.array\ |
|
| 192 |
# ([[float(f) for f in line.split()[0:3]]])) |
|
| 193 |
tmp=np.array([[float(f) for f in line.split()[0:3]]]) |
|
| 194 |
forces=np.concatenate((forces,tmp)) |
|
| 195 |
#note the '-' sign for turbomole, to get forces |
|
| 196 |
self.forces = (-np.delete(forces, np.s_[0:1], axis=0))*Hartree/Bohr |
|
Formats disponibles : Unified diff