Coverage for /builds/debichem-team/python-ase/ase/calculators/counterions.py : 96.72%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import numpy as np
2from ase.calculators.calculator import Calculator
3from ase import units
5k_c = units.Hartree * units.Bohr
8class AtomicCounterIon(Calculator):
9 implemented_properties = ['energy', 'forces']
11 def __init__(self, charge, epsilon, sigma, sites_per_mol=1,
12 rc=7.0, width=1.0):
13 """ Counter Ion Calculator.
15 A very simple, nonbonded (Coulumb and LJ)
16 interaction calculator meant for single atom ions
17 to charge neutralize systems (and nothing else)...
18 """
19 self.rc = rc
20 self.width = width
21 self.sites_per_mol = sites_per_mol
22 self.epsilon = epsilon
23 self.sigma = sigma
24 self.charge = charge
25 Calculator.__init__(self)
27 def add_virtual_sites(self, positions):
28 return positions
30 def get_virtual_charges(self, atoms):
31 charges = np.tile(self.charge, len(atoms) // self.sites_per_mol)
32 return charges
34 def redistribute_forces(self, forces):
35 return forces
37 def calculate(self, atoms, properties, system_changes):
38 Calculator.calculate(self, atoms, properties, system_changes)
40 R = atoms.get_positions()
41 charges = self.get_virtual_charges(atoms)
42 pbc = atoms.pbc
44 energy = 0.0
45 forces = np.zeros_like(atoms.get_positions())
47 for m in range(len(atoms)):
48 D = R[m + 1:] - R[m]
49 shift = np.zeros_like(D)
50 for i, periodic in enumerate(pbc):
51 if periodic:
52 L = atoms.cell.diagonal()[i]
53 shift[:, i] = (D[:, i] + L / 2) % L - L / 2 - D[:, i]
54 D += shift
55 d2 = (D**2).sum(1)
56 d = d2**0.5
58 x1 = d > self.rc - self.width
59 x2 = d < self.rc
60 x12 = np.logical_and(x1, x2)
61 y = (d[x12] - self.rc + self.width) / self.width
62 t = np.zeros(len(d)) # cutoff function
63 t[x2] = 1.0
64 t[x12] -= y**2 * (3.0 - 2.0 * y)
65 dtdd = np.zeros(len(d))
66 dtdd[x12] -= 6.0 / self.width * y * (1.0 - y)
68 c6 = (self.sigma**2 / d2)**3
69 c12 = c6**2
70 e_lj = 4 * self.epsilon * (c12 - c6)
71 e_c = k_c * charges[m + 1:] * charges[m] / d
73 energy += np.dot(t, e_lj)
74 energy += np.dot(t, e_c)
76 F = (24 * self.epsilon * (2 * c12 - c6) / d2 * t -
77 e_lj * dtdd / d)[:, None] * D
79 forces[m] -= F.sum(0)
80 forces[m + 1:] += F
82 F = (e_c / d2 * t)[:, None] * D \
83 - (e_c * dtdd / d)[:, None] * D
85 forces[m] -= F.sum(0)
86 forces[m + 1:] += F
88 self.results['energy'] = energy
89 self.results['forces'] = forces