From 48a2e58989e3f7dabc9c4a45dd1120957b102779 Mon Sep 17 00:00:00 2001 From: Mariana Rossi 2 Date: Sun, 10 Sep 2023 21:10:23 +0200 Subject: [PATCH] This commit only "blacks" the changes. --- ipi/engine/forcefields.py | 322 ++++++++++++++++++++------------------ 1 file changed, 171 insertions(+), 151 deletions(-) diff --git a/ipi/engine/forcefields.py b/ipi/engine/forcefields.py index 574f076c9..8e53817eb 100644 --- a/ipi/engine/forcefields.py +++ b/ipi/engine/forcefields.py @@ -81,13 +81,13 @@ class ForceField(dobject): """ def __init__( - self, - latency=1.0, - name="", - pars=None, - dopbc=True, - active=np.array([-1]), - threaded=False, + self, + latency=1.0, + name="", + pars=None, + dopbc=True, + active=np.array([-1]), + threaded=False, ): """Initialises ForceField. @@ -315,14 +315,14 @@ class FFSocket(ForceField): """ def __init__( - self, - latency=1.0, - name="", - pars=None, - dopbc=True, - active=np.array([-1]), - threaded=True, - interface=None, + self, + latency=1.0, + name="", + pars=None, + dopbc=True, + active=np.array([-1]), + threaded=True, + interface=None, ): """Initialises FFSocket. @@ -428,10 +428,10 @@ def evaluate(self, r): f = np.zeros(q.shape) for i in range(1, nat): dij = q[i] - q[:i] - rij2 = (dij ** 2).sum(axis=1) + rij2 = (dij**2).sum(axis=1) x6 = (self.sigma2 / rij2) ** 3 - x12 = x6 ** 2 + x12 = x6**2 v += (x12 - x6).sum() dij *= (self.sixepsfour * (2.0 * x12 - x6) / rij2)[:, np.newaxis] @@ -461,16 +461,16 @@ class FFdmd(ForceField): """ def __init__( - self, - latency=1.0e-3, - name="", - coupling=None, - freq=0.0, - dtdmd=0.0, - dmdstep=0, - pars=None, - dopbc=False, - threaded=False, + self, + latency=1.0e-3, + name="", + coupling=None, + freq=0.0, + dtdmd=0.0, + dmdstep=0, + pars=None, + dopbc=False, + threaded=False, ): """Initialises FFdmd. @@ -532,7 +532,7 @@ def evaluate(self, r): # rij = np.sqrt((dij ** 2).sum(axis=1)) # KF's implementation: dij, rij = vector_separation(cell_h, cell_ih, q[i], q[:i]) - cij = self.coupling[i * (i - 1) // 2: i * (i + 1) // 2] + cij = self.coupling[i * (i - 1) // 2 : i * (i + 1) // 2] prefac = np.dot( cij, rij ) # for each i it has the distances to all indexes previous @@ -574,15 +574,15 @@ class FFDebye(ForceField): """ def __init__( - self, - latency=1.0, - name="", - H=None, - xref=None, - vref=0.0, - pars=None, - dopbc=False, - threaded=False, + self, + latency=1.0, + name="", + H=None, + xref=None, + vref=0.0, + pars=None, + dopbc=False, + threaded=False, ): """Initialises FFDebye. @@ -661,15 +661,15 @@ class FFPlumed(ForceField): """ def __init__( - self, - latency=1.0e-3, - name="", - pars=None, - dopbc=False, - threaded=False, - init_file="", - plumeddat="", - plumedstep=0, + self, + latency=1.0e-3, + name="", + pars=None, + dopbc=False, + threaded=False, + init_file="", + plumeddat="", + plumedstep=0, ): """Initialises FFPlumed. @@ -817,21 +817,21 @@ class FFYaff(ForceField): """Use Yaff as a library to construct a force field""" def __init__( - self, - latency=1.0, - name="", - threaded=False, - yaffpara=None, - yaffsys=None, - yafflog="yaff.log", - rcut=18.89726133921252, - alpha_scale=3.5, - gcut_scale=1.1, - skin=0, - smooth_ei=False, - reci_ei="ewald", - pars=None, - dopbc=False, + self, + latency=1.0, + name="", + threaded=False, + yaffpara=None, + yaffsys=None, + yafflog="yaff.log", + rcut=18.89726133921252, + alpha_scale=3.5, + gcut_scale=1.1, + skin=0, + smooth_ei=False, + reci_ei="ewald", + pars=None, + dopbc=False, ): """Initialises FFYaff and enables a basic Yaff force field. @@ -935,13 +935,13 @@ class FFsGDML(ForceField): """ def __init__( - self, - latency=1.0, - name="", - threaded=False, - sGDML_model=None, - pars=None, - dopbc=False, + self, + latency=1.0, + name="", + threaded=False, + sGDML_model=None, + pars=None, + dopbc=False, ): """Initialises FFsGDML @@ -1048,21 +1048,21 @@ class FFCommittee(ForceField): committee of potentials, and implements the weighted baseline method.""" def __init__( - self, - latency=1.0, - name="", - pars=None, - dopbc=True, - active=np.array([-1]), - threaded=True, - fflist=[], - ffweights=[], - alpha=1.0, - baseline_name="", - baseline_uncertainty=-1.0, - active_thresh=0.0, - active_out=None, - comm_type="default", + self, + latency=1.0, + name="", + pars=None, + dopbc=True, + active=np.array([-1]), + threaded=True, + fflist=[], + ffweights=[], + alpha=1.0, + baseline_name="", + baseline_uncertainty=-1.0, + active_thresh=0.0, + active_out=None, + comm_type="default", ): # force threaded mode as otherwise it cannot have threaded children super(FFCommittee, self).__init__( @@ -1139,8 +1139,13 @@ def check_finish(self, r): def gather(self, r): """Collects results from all sub-requests, and assemble the committee of models.""" - if (self.comm_type == "default"): - r["result"] = [0.0, np.zeros(len(r["pos"]), float), np.zeros((3, 3), float), ""] + if self.comm_type == "default": + r["result"] = [ + 0.0, + np.zeros(len(r["pos"]), float), + np.zeros((3, 3), float), + "", + ] # list of pointers to the forcefield requests. shallow copy so we can remove stuff com_handles = r["ff_handles"].copy() @@ -1164,14 +1169,19 @@ def gather(self, r): xtrs = [ff_r["result"][3] for ff_r in com_handles] print("quants normal ", pots, frcs, virs) - elif (self.comm_type == "single-extras"): + elif self.comm_type == "single-extras": # Check that we indeed have a single ff - if (len(self.fflist) != 1 and self.baseline_name == ""): + if len(self.fflist) != 1 and self.baseline_name == "": raise ValueError( "This type of committee expects a single socket, or one socket and one baseline" ) - r["result"] = [0.0, np.zeros(len(r["pos"]), float), np.zeros((3, 3), float), ""] + r["result"] = [ + 0.0, + np.zeros(len(r["pos"]), float), + np.zeros((3, 3), float), + "", + ] # list of pointers to the forcefield requests. shallow copy so we can remove stuff com_handles = r["ff_handles"].copy() @@ -1192,10 +1202,20 @@ def gather(self, r): # debug for ff_r in com_handles: print("dict", ff_r["result"][3]) - pots = np.squeeze(np.array([ff_r["result"][3]["committee_pot"] for ff_r in com_handles])) - frcs = np.squeeze(np.array([ff_r["result"][3]["committee_force"] for ff_r in com_handles])) - virs = np.squeeze(np.array([ff_r["result"][3]["committee_virial"] for ff_r in com_handles])) - xtrs = [ff_r["result"][3] for ff_r in com_handles] # MR - no extras here kinda. not sure yet about this + pots = np.squeeze( + np.array([ff_r["result"][3]["committee_pot"] for ff_r in com_handles]) + ) + frcs = np.squeeze( + np.array([ff_r["result"][3]["committee_force"] for ff_r in com_handles]) + ) + virs = np.squeeze( + np.array( + [ff_r["result"][3]["committee_virial"] for ff_r in com_handles] + ) + ) + xtrs = [ + ff_r["result"][3] for ff_r in com_handles + ] # MR - no extras here kinda. not sure yet about this # debug # print("quants0 ", pots, frcs, virs, pots.shape, frcs.shape, virs.shape, pots.dtype, frcs.dtype, virs.dtype) # MR: from now on everything else should be the same, hopefully @@ -1233,43 +1253,43 @@ def gather(self, r): # and V_committee the committee error. Then # V = V_baseline + s_b^2/(s_c^2+s_b^2) V_committe - s_b2 = self.baseline_uncertainty ** 2 + s_b2 = self.baseline_uncertainty**2 nmodels = len(pots) uncertain_frc = ( - self.alpha ** 2 - * np.sum( - [ - (pot - mean_pot) * (frc - mean_frc) - for pot, frc in zip(pots, frcs) - ], - axis=0, - ) - / (nmodels - 1) + self.alpha**2 + * np.sum( + [ + (pot - mean_pot) * (frc - mean_frc) + for pot, frc in zip(pots, frcs) + ], + axis=0, + ) + / (nmodels - 1) ) uncertain_vir = ( - self.alpha ** 2 - * np.sum( - [ - (pot - mean_pot) * (vir - mean_vir) - for pot, vir in zip(pots, virs) - ], - axis=0, - ) - / (nmodels - 1) + self.alpha**2 + * np.sum( + [ + (pot - mean_pot) * (vir - mean_vir) + for pot, vir in zip(pots, virs) + ], + axis=0, + ) + / (nmodels - 1) ) # Computes the final average energetics final_pot = baseline_pot + mean_pot * s_b2 / (s_b2 + var_pot) final_frc = ( - baseline_frc - + mean_frc * s_b2 / (s_b2 + var_pot) - - 2.0 * mean_pot * s_b2 / (s_b2 + var_pot) ** 2 * uncertain_frc + baseline_frc + + mean_frc * s_b2 / (s_b2 + var_pot) + - 2.0 * mean_pot * s_b2 / (s_b2 + var_pot) ** 2 * uncertain_frc ) final_vir = ( - baseline_vir - + mean_vir * s_b2 / (s_b2 + var_pot) - - 2.0 * mean_pot * s_b2 / (s_b2 + var_pot) ** 2 * uncertain_vir + baseline_vir + + mean_vir * s_b2 / (s_b2 + var_pot) + - 2.0 * mean_pot * s_b2 / (s_b2 + var_pot) ** 2 * uncertain_vir ) # Sets the output of the committee model. @@ -1392,10 +1412,10 @@ def init_photon(self): # construct varepsilon array for all photons self.varepsilon_k = self.E0 self.varepsilon_klambda = ( - self.E0 * self.omega_klambda / np.min(self.omega_klambda) + self.E0 * self.omega_klambda / np.min(self.omega_klambda) ) self.varepsilon_klambda3 = ( - self.E0 * self.omega_klambda3 / np.min(self.omega_klambda3) + self.E0 * self.omega_klambda3 / np.min(self.omega_klambda3) ) # cavity mode function acting on the dipole moment @@ -1416,7 +1436,7 @@ def split_atom_ph_coord(self, pos): """ if self.apply_photon: pos_at = pos[: -self.n_photon_3] - pos_ph = pos[-self.n_photon_3:] + pos_ph = pos[-self.n_photon_3 :] self.pos_ph = pos_ph else: pos_at = pos @@ -1436,7 +1456,7 @@ def get_ph_energy(self, dx_array, dy_array): total energy of photonic system """ # calculate the photonic potential energy - e_ph = np.sum(0.5 * self.omega_klambda3 ** 2 * self.pos_ph ** 2) + e_ph = np.sum(0.5 * self.omega_klambda3**2 * self.pos_ph**2) # calculate the dot products between mode functions and dipole array d_dot_f_x = np.dot(self.ftilde_kx, dx_array) @@ -1445,10 +1465,10 @@ def get_ph_energy(self, dx_array, dy_array): # calculate the light-matter interaction if self.ph_rep == "loose": e_int_x = np.sum( - self.varepsilon_k * d_dot_f_x * self.pos_ph[: self.n_mode * 3: 3] + self.varepsilon_k * d_dot_f_x * self.pos_ph[: self.n_mode * 3 : 3] ) e_int_y = np.sum( - self.varepsilon_k * d_dot_f_y * self.pos_ph[1 + self.n_mode * 3:: 3] + self.varepsilon_k * d_dot_f_y * self.pos_ph[1 + self.n_mode * 3 :: 3] ) elif self.ph_rep == "dense": e_int_x = np.sum(self.varepsilon_k * d_dot_f_x * self.pos_ph[::3]) @@ -1456,8 +1476,8 @@ def get_ph_energy(self, dx_array, dy_array): # calculate the dipole self-energy term dse = np.sum( - (self.varepsilon_k ** 2 / 2.0 / self.omega_k ** 2) - * (d_dot_f_x ** 2 + d_dot_f_y ** 2) + (self.varepsilon_k**2 / 2.0 / self.omega_k**2) + * (d_dot_f_x**2 + d_dot_f_y**2) ) e_tot = e_ph + e_int_x + e_int_y + dse @@ -1476,7 +1496,7 @@ def get_ph_forces(self, dx_array, dy_array): force array of all photonic dimensions (3*nphoton) [1x, 1y, 1z, 2x..] """ # calculat the bare photonic contribution of the force - f_ph = -self.omega_klambda3 ** 2 * self.pos_ph + f_ph = -self.omega_klambda3**2 * self.pos_ph # calculate the dot products between mode functions and dipole array d_dot_f_x = np.dot(self.ftilde_kx, dx_array) @@ -1484,8 +1504,8 @@ def get_ph_forces(self, dx_array, dy_array): # calculate the force due to light-matter interactions if self.ph_rep == "loose": - f_ph[: self.n_mode * 3: 3] -= self.varepsilon_k * d_dot_f_x - f_ph[self.n_mode * 3 + 1:: 3] -= self.varepsilon_k * d_dot_f_y + f_ph[: self.n_mode * 3 : 3] -= self.varepsilon_k * d_dot_f_x + f_ph[self.n_mode * 3 + 1 :: 3] -= self.varepsilon_k * d_dot_f_y elif self.ph_rep == "dense": f_ph[::3] -= self.varepsilon_k * d_dot_f_x f_ph[1::3] -= self.varepsilon_k * d_dot_f_y @@ -1510,13 +1530,13 @@ def get_nuc_cav_forces(self, dx_array, dy_array, charge_array_bath): # cavity force on x direction if self.ph_rep == "loose": - Ekx = self.varepsilon_k * self.pos_ph[: self.n_mode * 3: 3] - Eky = self.varepsilon_k * self.pos_ph[self.n_mode * 3 + 1:: 3] + Ekx = self.varepsilon_k * self.pos_ph[: self.n_mode * 3 : 3] + Eky = self.varepsilon_k * self.pos_ph[self.n_mode * 3 + 1 :: 3] elif self.ph_rep == "dense": Ekx = self.varepsilon_k * self.pos_ph[::3] Eky = self.varepsilon_k * self.pos_ph[1::3] - Ekx += self.varepsilon_k ** 2 / self.omega_k ** 2 * d_dot_f_x - Eky += self.varepsilon_k ** 2 / self.omega_k ** 2 * d_dot_f_y + Ekx += self.varepsilon_k**2 / self.omega_k**2 * d_dot_f_x + Eky += self.varepsilon_k**2 / self.omega_k**2 * d_dot_f_y # dimension of independent baths (xy grid points) coeff_x = np.dot(np.transpose(Ekx), self.ftilde_kx) @@ -1536,19 +1556,19 @@ class FFCavPhSocket(FFSocket): """ def __init__( - self, - latency=1.0, - name="", - pars=None, - dopbc=False, - active=np.array([-1]), - threaded=True, - interface=None, - charge_array=None, - apply_photon=True, - E0=1e-4, - omega_c=0.01, - ph_rep="loose", + self, + latency=1.0, + name="", + pars=None, + dopbc=False, + active=np.array([-1]), + threaded=True, + interface=None, + charge_array=None, + apply_photon=True, + E0=1e-4, + omega_c=0.01, + ph_rep="loose", ): """Initialises FFCavPhFPSocket. @@ -1613,7 +1633,7 @@ def calc_dipole_xyz_mm(self, pos, n_bath, charge_array_bath): dx_array, dy_array, dz_array = [], [], [] for idx in range(n_bath): - pos_bath = pos[ndim_local * idx: ndim_local * (idx + 1)] + pos_bath = pos[ndim_local * idx : ndim_local * (idx + 1)] # check the dimension of charge array if np.size(pos_bath[::3]) != np.size(charge_array_bath): softexit.trigger( @@ -1696,8 +1716,8 @@ def queue(self, atoms, cell, reqid=-1): # 2. for atomic coordinates, we now evaluate their atomic forces for idx in range(self.n_independent_bath): pbcpos_local = pbcpos_atoms[ - ndim_local * idx: ndim_local * (idx + 1) - ].copy() + ndim_local * idx : ndim_local * (idx + 1) + ].copy() iactive_local = self.iactive[0:ndim_local] # Let's try to do PBC for the small regions if self.dopbc: @@ -1763,7 +1783,7 @@ def queue(self, atoms, cell, reqid=-1): for idx, newreq in enumerate(newreq_lst): u, f, vir, extra = newreq["result"] result_tot[0] += u - result_tot[1][ndim_local * idx: ndim_local * (idx + 1)] = f + result_tot[1][ndim_local * idx : ndim_local * (idx + 1)] = f result_tot[2] += vir result_tot[3][idx] = extra @@ -1776,8 +1796,8 @@ def queue(self, atoms, cell, reqid=-1): ) # check the size of photon modes + molecules to match the total number of particles if ( - self.ph.n_photon + self.n_independent_bath * self.charge_array.size - != int(len(pbcpos) // 3) + self.ph.n_photon + self.n_independent_bath * self.charge_array.size + != int(len(pbcpos) // 3) ): softexit.trigger( "Total number of photons + molecules does not match total number of particles"