Révision 3 ase/calculators/qmx.py
| qmx.py (revision 3) | ||
|---|---|---|
| 8 | 8 |
import numpy as np |
| 9 | 9 |
from general import Calculator |
| 10 | 10 |
from ase.embed import Embed |
| 11 |
from ase.units import Hartree |
|
| 11 | 12 |
|
| 13 |
from copy import deepcopy |
|
| 14 |
|
|
| 12 | 15 |
import sys, os |
| 13 | 16 |
|
| 14 | 17 |
class Qmx(Calculator): |
| 15 |
def __init__(self, calculator_low, calculator_high): |
|
| 16 |
self.string_params = {}
|
|
| 17 |
|
|
| 18 |
self.forces=np.zeros((2,2)) |
|
| 18 |
def __init__(self, calculator_high, calculator_low): |
|
| 19 | 19 |
self._constraints=None |
| 20 |
|
|
| 21 |
self.calculator_low = calculator_low
|
|
| 22 |
self.calculator_high = calculator_high |
|
| 23 |
|
|
| 20 |
self.calculator_low_cluster = calculator_low |
|
| 21 |
self.calculator_low_system = deepcopy(calculator_low)
|
|
| 22 |
self.calculator_high = calculator_high
|
|
| 23 |
|
|
| 24 | 24 |
def get_energy_subsystem(self, path, calculator, atoms, force_consistent): |
| 25 |
# go to directory and calculate energies |
|
| 26 |
print "running energy in: ", path |
|
| 25 | 27 |
os.chdir(path) |
| 26 |
calculator.set_atoms(atoms)
|
|
| 27 |
energy = calculator.get_potential_energy(atoms)
|
|
| 28 |
atoms.set_calculator(calculator)
|
|
| 29 |
energy = atoms.get_potential_energy()
|
|
| 28 | 30 |
os.chdir("..")
|
| 29 | 31 |
return energy |
| 30 |
|
|
| 32 |
|
|
| 31 | 33 |
def get_forces_subsystem(self, path, calculator, atoms): |
| 34 |
# go to directory and calculate forces |
|
| 35 |
print "running forces in: ", path |
|
| 32 | 36 |
os.chdir(path) |
| 33 |
calculator.set_atoms(atoms)
|
|
| 34 |
forces = calculator.get_forces(atoms)
|
|
| 37 |
atoms.set_calculator(calculator)
|
|
| 38 |
forces = atoms.get_forces()
|
|
| 35 | 39 |
os.chdir("..")
|
| 36 | 40 |
return forces |
| 37 |
|
|
| 38 |
|
|
| 41 |
|
|
| 39 | 42 |
def get_potential_energy(self, embed, force_consistent=False): |
| 40 |
# perform energy calculations
|
|
| 41 |
e_sys_lo = self.get_energy_subsystem("system.low-level", self.calculator_low, embed.get_system(), force_consistent)
|
|
| 42 |
e_cl_lo = self.get_energy_subsystem("cluster.low-level", self.calculator_low, embed.get_cluster(), force_consistent)
|
|
| 43 |
# perform energy calculations |
|
| 44 |
e_sys_lo = self.get_energy_subsystem("system.low-level", self.calculator_low_system, embed.get_system(), force_consistent)
|
|
| 45 |
e_cl_lo = self.get_energy_subsystem("cluster.low-level", self.calculator_low_cluster, embed.get_cluster(), force_consistent)
|
|
| 43 | 46 |
e_cl_hi = self.get_energy_subsystem("cluster.high-level", self.calculator_high, embed.get_cluster(), force_consistent)
|
| 44 | 47 |
# calculate energies |
| 45 | 48 |
energy = e_sys_lo - e_cl_lo + e_cl_hi |
| 49 |
# print energies |
|
| 46 | 50 |
print "%20s=%15s - %15s + %15s" %("E(C:S)", "E(S-MM)", "E(C-MM)", "E(C-QM)")
|
| 47 | 51 |
print "%20f=%15f - %15f + %15f" %(energy, e_sys_lo, e_cl_lo, e_cl_hi) |
| 52 |
# set energies and return |
|
| 48 | 53 |
if force_consistent: |
| 49 | 54 |
self.energy_free = energy |
| 50 | 55 |
return self.energy_free |
| ... | ... | |
| 55 | 60 |
def get_forces(self, embed): |
| 56 | 61 |
atom_map_sys_cl = embed.atom_map_sys_cl |
| 57 | 62 |
# get forces for the three systems |
| 58 |
f_sys_lo = self.get_forces_subsystem("system.low-level", self.calculator_low, embed.get_system())
|
|
| 59 |
f_cl_lo = self.get_forces_subsystem("cluster.low-level", self.calculator_low, embed.get_cluster())
|
|
| 63 |
f_sys_lo = self.get_forces_subsystem("system.low-level", self.calculator_low_system, embed.get_system())
|
|
| 64 |
f_cl_lo = self.get_forces_subsystem("cluster.low-level", self.calculator_low_cluster, embed.get_cluster())
|
|
| 60 | 65 |
f_cl_hi = self.get_forces_subsystem("cluster.high-level", self.calculator_high, embed.get_cluster())
|
| 61 | 66 |
# forces correction for the atoms |
| 62 | 67 |
f_cl = f_cl_hi - f_cl_lo |
| 63 |
#number of atoms |
|
| 68 |
# number of atoms
|
|
| 64 | 69 |
nat_sys = len(embed) |
| 65 | 70 |
# lo-sys + (hi-lo) |
| 66 | 71 |
for iat_sys in xrange(nat_sys): |
| 67 | 72 |
iat_cl = atom_map_sys_cl[iat_sys] |
| 68 | 73 |
if iat_cl > -1: |
| 69 | 74 |
f_sys_lo[iat_sys] += f_cl[iat_cl] |
| 70 |
|
|
| 75 |
|
|
| 71 | 76 |
# correct gradients |
| 72 | 77 |
# Reference: Eichler, Koelmel, Sauer, J. of Comput. Chem., 18(4). 1997, 463-477. |
| 73 |
for cell_L, iat_cl, iat_sys, r, iat_link in embed.linkatoms: |
|
| 78 |
for cell_L, iat_cl_sys, iat_sys, r, iat_link in embed.linkatoms:
|
|
| 74 | 79 |
# calculate the bond distance (r_bond) at the border |
| 75 |
xyz = embed[iat_sys].get_position() - embed.get_cluster()[iat_cl].get_position() + cell_L
|
|
| 80 |
xyz = embed[iat_sys].get_position() - embed[iat_cl_sys].get_position() + cell_L
|
|
| 76 | 81 |
|
| 77 | 82 |
# calculate the bond lenght and the factor f |
| 78 | 83 |
rbond = np.sqrt(np.dot(xyz, xyz)) |
| 79 | 84 |
f = r / rbond |
| 80 | 85 |
#normalize xyz |
| 81 | 86 |
xyz /= rbond |
| 82 |
|
|
| 83 | 87 |
# receive the gradients for the link atom |
| 84 |
fH = f_cl[iat_link]
|
|
| 85 |
# Skalarprodukt fH, xyz
|
|
| 86 |
fs = np.dot(xyz, fH)
|
|
| 87 |
|
|
| 88 |
fL = f_cl[iat_link]
|
|
| 89 |
# dot product fL, xyz
|
|
| 90 |
fs = np.dot(xyz, fL)
|
|
| 91 |
# apply corrections for each direction |
|
| 88 | 92 |
for idir in xrange(3): |
| 89 | 93 |
# correct the atom in the system |
| 90 |
f_sys_lo[iat_sys][idir] += f*fH[idir] - f*fs*xyz[idir]
|
|
| 94 |
f_sys_lo[iat_sys][idir] += f*fL[idir] - f*fs*xyz[idir]
|
|
| 91 | 95 |
# correct the atom in the cluster |
| 92 |
f_sys_lo[iat_cl][idir] += (1-f)*fH[idir] + f*fs*xyz[idir]
|
|
| 96 |
f_sys_lo[iat_cl_sys][idir] += (1-f)*fL[idir] + f*fs*xyz[idir]
|
|
| 93 | 97 |
return f_sys_lo |
| 94 |
|
|
| 95 |
def set_atoms(self, atoms): |
|
| 96 |
return |
|
| 97 |
|
|
| 98 |
|
|
| 99 |
|
|
Formats disponibles : Unified diff