Révision 5 ase/calculators/qmx.py
| qmx.py (revision 5) | ||
|---|---|---|
| 15 | 15 |
import sys, os |
| 16 | 16 |
|
| 17 | 17 |
class Qmx(Calculator): |
| 18 |
def __init__(self, calculator_high, calculator_low):
|
|
| 18 |
def __init__(self, calculator_high_cluster, calculator_low_cluster, calculator_low_system=None, print_forces=False):
|
|
| 19 | 19 |
self._constraints=None |
| 20 |
self.calculator_low_cluster = calculator_low |
|
| 21 |
self.calculator_low_system = deepcopy(calculator_low) |
|
| 22 |
self.calculator_high = calculator_high |
|
| 20 |
|
|
| 21 |
self.calculator_low_cluster = calculator_low_cluster |
|
| 22 |
self.calculator_low_system = calculator_low_system |
|
| 23 |
if self.calculator_low_system is None: |
|
| 24 |
self.calculator_low_system = deepcopy(calculator_low_cluster) |
|
| 25 |
self.calculator_high_cluster = calculator_high_cluster |
|
| 26 |
self.print_forces = print_forces |
|
| 23 | 27 |
|
| 24 | 28 |
def get_energy_subsystem(self, path, calculator, atoms, force_consistent): |
| 25 | 29 |
# go to directory and calculate energies |
| ... | ... | |
| 43 | 47 |
# perform energy calculations |
| 44 | 48 |
e_sys_lo = self.get_energy_subsystem("system.low-level", self.calculator_low_system, embed.get_system(), force_consistent)
|
| 45 | 49 |
e_cl_lo = self.get_energy_subsystem("cluster.low-level", self.calculator_low_cluster, embed.get_cluster(), force_consistent)
|
| 46 |
e_cl_hi = self.get_energy_subsystem("cluster.high-level", self.calculator_high, embed.get_cluster(), force_consistent)
|
|
| 50 |
e_cl_hi = self.get_energy_subsystem("cluster.high-level", self.calculator_high_cluster, embed.get_cluster(), force_consistent)
|
|
| 47 | 51 |
# calculate energies |
| 48 | 52 |
energy = e_sys_lo - e_cl_lo + e_cl_hi |
| 49 | 53 |
# print energies |
| 50 |
print "%20s=%15s - %15s + %15s" %("E(C:S)", "E(S-MM)", "E(C-MM)", "E(C-QM)")
|
|
| 51 |
print "%20f=%15f - %15f + %15f" %(energy, e_sys_lo, e_cl_lo, e_cl_hi)
|
|
| 54 |
print "%20s = %15s - %15s + %15s" %("E(C:S)", "E(S-MM)", "E(C-MM)", "E(C-QM)")
|
|
| 55 |
print "%20f = %15f - %15f + %15f" %(energy, e_sys_lo, e_cl_lo, e_cl_hi)
|
|
| 52 | 56 |
# set energies and return |
| 53 | 57 |
if force_consistent: |
| 54 | 58 |
self.energy_free = energy |
| ... | ... | |
| 62 | 66 |
# get forces for the three systems |
| 63 | 67 |
f_sys_lo = self.get_forces_subsystem("system.low-level", self.calculator_low_system, embed.get_system())
|
| 64 | 68 |
f_cl_lo = self.get_forces_subsystem("cluster.low-level", self.calculator_low_cluster, embed.get_cluster())
|
| 65 |
f_cl_hi = self.get_forces_subsystem("cluster.high-level", self.calculator_high, embed.get_cluster())
|
|
| 69 |
f_cl_hi = self.get_forces_subsystem("cluster.high-level", self.calculator_high_cluster, embed.get_cluster())
|
|
| 70 |
|
|
| 66 | 71 |
# forces correction for the atoms |
| 67 | 72 |
f_cl = f_cl_hi - f_cl_lo |
| 68 |
# number of atoms |
|
| 69 |
nat_sys = len(embed) |
|
| 73 |
|
|
| 74 |
if self.print_forces: |
|
| 75 |
cluster=embed.get_cluster() |
|
| 76 |
print "Forces: System LOW - Cluster LOW + Cluster HIGH" |
|
| 77 |
for iat_sys in xrange(len(embed)): |
|
| 78 |
print "%-2s (" % embed[iat_sys].get_symbol(),
|
|
| 79 |
for idir in xrange(3): |
|
| 80 |
print "%10.6f" % f_sys_lo[iat_sys][idir], |
|
| 81 |
print ") <system LOW>" |
|
| 82 |
|
|
| 83 |
iat_cl = atom_map_sys_cl[iat_sys] |
|
| 84 |
if iat_cl > -1: |
|
| 85 |
print "%s" % "- (",
|
|
| 86 |
for idir in xrange(3): |
|
| 87 |
print "%10.6f" % f_cl_lo[iat_cl][idir], |
|
| 88 |
print ") <cluster LOW>" |
|
| 89 |
print "%s" % "+ (",
|
|
| 90 |
for idir in xrange(3): |
|
| 91 |
print "%10.6f" % f_cl_hi[iat_cl][idir], |
|
| 92 |
print ") <cluster HIGH>" |
|
| 93 |
|
|
| 94 |
|
|
| 95 |
|
|
| 70 | 96 |
# lo-sys + (hi-lo) |
| 71 |
for iat_sys in xrange(nat_sys):
|
|
| 97 |
for iat_sys in xrange(len(embed)):
|
|
| 72 | 98 |
iat_cl = atom_map_sys_cl[iat_sys] |
| 73 | 99 |
if iat_cl > -1: |
| 74 | 100 |
f_sys_lo[iat_sys] += f_cl[iat_cl] |
| 75 |
|
|
| 101 |
# some settings for the output |
|
| 102 |
i_change = np.zeros(len(embed), int) |
|
| 103 |
if self.print_forces: |
|
| 104 |
f_sys_lo_orig = f_sys_lo.copy() |
|
| 76 | 105 |
# correct gradients |
| 77 | 106 |
# Reference: Eichler, Koelmel, Sauer, J. of Comput. Chem., 18(4). 1997, 463-477. |
| 78 | 107 |
for cell_L, iat_cl_sys, iat_sys, r, iat_link in embed.linkatoms: |
| 79 | 108 |
# calculate the bond distance (r_bond) at the border |
| 80 | 109 |
xyz = embed[iat_sys].get_position() - embed[iat_cl_sys].get_position() + cell_L |
| 81 |
|
|
| 82 | 110 |
# calculate the bond lenght and the factor f |
| 83 | 111 |
rbond = np.sqrt(np.dot(xyz, xyz)) |
| 84 | 112 |
f = r / rbond |
| ... | ... | |
| 87 | 115 |
# receive the gradients for the link atom |
| 88 | 116 |
fL = f_cl[iat_link] |
| 89 | 117 |
# dot product fL, xyz |
| 90 |
fs = np.dot(xyz, fL)
|
|
| 118 |
fs = np.dot(fL, xyz)
|
|
| 91 | 119 |
# apply corrections for each direction |
| 120 |
i_change[iat_sys] = 1 |
|
| 121 |
i_change[iat_cl_sys] = 1 |
|
| 92 | 122 |
for idir in xrange(3): |
| 93 | 123 |
# correct the atom in the system |
| 94 |
f_sys_lo[iat_sys][idir] += f*fL[idir] - f*fs*xyz[idir]
|
|
| 124 |
f_sys_lo[iat_sys][idir] = f_sys_lo[iat_sys][idir] + f*fL[idir] - f*fs*xyz[idir]
|
|
| 95 | 125 |
# correct the atom in the cluster |
| 96 |
f_sys_lo[iat_cl_sys][idir] += (1-f)*fL[idir] + f*fs*xyz[idir] |
|
| 126 |
f_sys_lo[iat_cl_sys][idir] = f_sys_lo[iat_cl_sys][idir] + (1-f)*fL[idir] + f*fs*xyz[idir] |
|
| 127 |
|
|
| 128 |
if self.print_forces: |
|
| 129 |
print " TOTAL FORCE (uncorrected : corrected) for link atoms" |
|
| 130 |
for iat_sys in xrange(len(embed)): |
|
| 131 |
print "%-2s (" % embed[iat_sys].get_symbol(),
|
|
| 132 |
for idir in xrange(3): |
|
| 133 |
print "%10.6f" % f_sys_lo_orig[iat_sys][idir], |
|
| 134 |
print ") : (",
|
|
| 135 |
for idir in xrange(3): |
|
| 136 |
print "%10.6f" % f_sys_lo[iat_sys][idir], |
|
| 137 |
print ")", |
|
| 138 |
if i_change[iat_sys]: |
|
| 139 |
print " *", |
|
| 140 |
|
|
| 141 |
|
|
| 142 |
|
|
| 97 | 143 |
return f_sys_lo |
Formats disponibles : Unified diff