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)

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