SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
COPSSImage.cpp
1 
21 #include "COPSSImage.h"
22 #include <math.h>
23 #include <algorithm>
24 #include <stdlib.h>
25 #include <stdio.h>
26 //#define DEBUG_SNAPSHOT_
27 //#define DEBUG_TMP_
28 namespace SSAGES
29 {
30 
32  {
33  //nothing to be done in presimulation
34  }
35 
37  {
38 
39  // Gather information
40  auto& types = snapshot->GetAtomTypes();
41  size_t nlocal_ = snapshot->GetNumAtoms();
42 
43  #ifdef DEBUG_SNAPSHOT_
44  fprintf(stdout, "Image method start runnning ... \n");
45 
46  auto& forces = snapshot->GetForces();
47  auto& positions = snapshot->GetPositions();
48  auto& charges = snapshot->GetCharges();
49 
50  fprintf(stdout, "check ids: ");
51  for (size_t i = 0; i < nlocal_; i++)
52  {
53  fprintf(stdout, "%d\n; ", ids[i]);
54  }
55 
56  fprintf(stdout, "check positions: ");
57  for (size_t i=0; i< nlocal_; i++)
58  {
59  for (size_t j=0; j<3; j++)
60  {
61  fprintf(stdout,"%f; ", positions[i][j]);
62  }
63  }
64 
65  fprintf(stdout, "\ncheck forces: ");
66  for (size_t i=0; i< nlocal_; i++)
67  {
68  for (size_t j=0; j<3; j++)
69  {
70  fprintf(stdout,"%f; ", forces[i][j]);
71  }
72  }
73 
74  fprintf(stdout, "\ncheck types: ");
75  for (size_t i=0; i< nlocal_; i++)
76  {
77  fprintf(stdout, "%d; ", types[i]);
78  }
79 
80  fprintf(stdout, "\ncheck charges: ");
81  for (size_t i=0; i< nlocal_; i++)
82  {
83  fprintf(stdout, "%f; ", charges[i]);
84  }
85 
86  fprintf(stdout, "\ncheck radius: ");
87  for (size_t i=0; i< nlocal_; i++)
88  {
89  fprintf(stdout,"%f; ", atomTypeRadius_[types[i]-1]);
90  }
91 
92  fprintf(stdout, "\ncheck nlocal: %d\n", nlocal_);
93  fprintf(stdout, "check eouter: %f\n", eouter_);
94  fprintf(stdout, "check qqrd2e: %f\n", qqrd2e_);
95  fprintf(stdout, "check einner: %f\n", einner_);
96  fprintf(stdout, "check ion-type-start: %d\n", ion_type_start_);
97  #endif
98 
99 
100  //iterate over all permutations {i;j;k}, j has to be polarizable
101 
102  for (size_t j=0; j < nlocal_; j++) if( types[j] < ion_type_start_ )
103  {
104  for(size_t i=0; i < nlocal_; i++) if(i != j)
105  {
106  for(size_t k=0; k < nlocal_; k++) if (k != j)
107  {
108  //Interaction: [I] <---- [J] <-------[k]
109  //F[i] = f{i;j;k}
110  //F[j] = -(f{i;j;k} + f{k;j;i})
111  //F[k] = f{k;j;i}
112  if (i == k)
113  {
114  force_pol(snapshot, i, j, k);
115  }
116  else if (i < k)
117  {
118  force_pol(snapshot, i, j, k);
119  force_pol(snapshot, k, j, i);
120  }
121  }
122  }
123  }
124 
125 
126  }
127 
129  {
130  //nothing to be done in post simulation
131  }
132 
133  // force acting on ith particle in a kernel {ith; jth; kth}
134  void COPSSImage::force_pol (Snapshot* snapshot,
135  size_t ith,
136  size_t jth,
137  size_t kth)
138  {
139  //fprintf(stdout, "compute: (%d, %d, %d)\n", ith, jth, kth);
140  auto& types_ = snapshot->GetAtomTypes();
141  auto& charges_ = snapshot->GetCharges();
142  auto& positions_ = snapshot->GetPositions();
143  auto& forces_ = snapshot->GetForces();
144  eouter_ = snapshot->GetDielectric();
145  qqrd2e_ = snapshot->Getqqrd2e();
146  aa = atomTypeRadius_[types_[jth]-1];
147  // helper variables
148  e_ = (1.0 - einner_ / eouter_) / (1.0 + einner_ / eouter_);
149  ginv_ = (1.0 + einner_ / eouter_);
150  //==================kernel function starts========================
151  // Rxkj, Rykj, Rzkj: vector points from j-th to k-th, in x, y, z direction respectively.
152  Rxkj = positions_[kth][0] - positions_[jth][0];
153  Rykj = positions_[kth][1] - positions_[jth][1];
154  Rzkj = positions_[kth][2] - positions_[jth][2];
155  // rkj: module of the vector between j-th and k-th
156  Rkj2 = Rxkj*Rxkj + Rykj*Rykj + Rzkj*Rzkj;
157  rkj = sqrt (Rkj2);
158  // ukj, vkj, wkj: unit vector points from j-th to k-th, in x, y, z direction repectively
159  ukj = Rxkj / rkj;
160  vkj = Rykj / rkj;
161  wkj = Rzkj / rkj;
162  // Rxij, Ryij, Rzij: vector points from i-th to j-th, in x, y, z direction respectively
163  Rxij = positions_[ith][0] - positions_[jth][0];
164  Ryij = positions_[ith][1] - positions_[jth][1];
165  Rzij = positions_[ith][2] - positions_[jth][2];
166  Rij2 = Rxij * Rxij + Ryij * Ryij + Rzij * Rzij;
167  rij = sqrt(Rij2);
168  // Auxiliary variables
169  aux1 = e_ * aa / rkj;
170  aux2 = aa * aa / rkj;
171  // Tmp variable for return value
172  std::vector<double> force_pol(3);
173  //===============delta_s,1 term of equation.29 in method paper
174  auxv_x_delta = Rxij - aux2 * ukj;
175  auxv_y_delta = Ryij - aux2 * vkj;
176  auxv_z_delta = Rzij - aux2 * wkj;
181  force_pol[0] = auxv_x_delta / aux3_delta;
182  force_pol[1] = auxv_y_delta / aux3_delta;
183  force_pol[2] = auxv_z_delta / aux3_delta;
184  //===============integration term of equation.29 in method paper
185  std::vector<double> xg_(5);
186  xlo = 0.0;
187  xhi = 1.0;
188  for (int ig = 0; ig < ngauss; ++ig)
189  {
190  xg_[ig] = 0.5 * (xhi - xlo) * xg0_[ig] + 0.5 * (xhi + xlo);
191  auxv_x_integ = Rxij - ( pow(xg_[ig] , ginv_) * aux2 * ukj );
192  auxv_y_integ = Ryij - ( pow(xg_[ig] , ginv_) * aux2 * vkj );
193  auxv_z_integ = Rzij - ( pow(xg_[ig] , ginv_) * aux2 * wkj );
195  auxv_y_integ * auxv_y_integ +
196  auxv_z_integ * auxv_z_integ);
197  aux3_integ = aux3Sqrt_integ * aux3Sqrt_integ * aux3Sqrt_integ;
198  force_pol[0] += 0.5 * (xhi-xlo) * (- auxv_x_integ / aux3_integ * wg0_[ig]);
199  force_pol[1] += 0.5 * (xhi-xlo) * (- auxv_y_integ / aux3_integ * wg0_[ig]);
200  force_pol[2] += 0.5 * (xhi-xlo) * (- auxv_z_integ / aux3_integ * wg0_[ig]);
201  }
202  // multiple by prefactor in eq.29
203  force_pol[0] *= ( 0.5 * qqrd2e_ * aux1 * charges_[ith] * charges_[kth] );
204  force_pol[1] *= ( 0.5 * qqrd2e_ * aux1 * charges_[ith] * charges_[kth] );
205  force_pol[2] *= ( 0.5 * qqrd2e_ * aux1 * charges_[ith] * charges_[kth] );
206 
207  for (size_t k = 0; k <3; ++k)
208  {
209  forces_[ith][k] += 2 * force_pol[k];
210  forces_[jth][k] -= 2 * force_pol[k];
211  }
212  } // end force_pol function
213 }//end namespace
214 
void PostIntegration(Snapshot *snapshot, const CVList &) override
Post-integration hook.
Definition: COPSSImage.cpp:36
double rij
Auxiliary variable r_ij.
Definition: COPSSImage.h:122
int ion_type_start_
Where non polarizable particles start.
Definition: COPSSImage.h:46
double Rykj
Auxiliary variable R_kj for image kernel functions (y component).
Definition: COPSSImage.h:98
double wg0_[5]
Magic numbers for the weight.
Definition: COPSSImage.h:67
double ukj
Auxiliary variable u_kj.
Definition: COPSSImage.h:125
unsigned GetNumAtoms() const
Get number of atoms in this snapshot.
Definition: Snapshot.h:200
double eouter_
Dielectric constant of outside continuum.
Definition: COPSSImage.h:49
std::vector< CollectiveVariable * > CVList
List of Collective Variables.
Definition: types.h:51
void PreSimulation(Snapshot *, const CVList &) override
Pre-simulation hook.
Definition: COPSSImage.cpp:31
double xhi
Upper value for x.
Definition: COPSSImage.h:58
const std::vector< Vector3 > & GetForces() const
Access the per-particle forces.
Definition: Snapshot.h:362
double Rxij
Auxiliary variable R_ij (x component).
Definition: COPSSImage.h:110
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:43
double aux3Sqrt_delta
Auxiliary variable, square root of delta.
Definition: COPSSImage.h:167
double auxv_x_integ
Auxiliary variable for integration (x component).
Definition: COPSSImage.h:140
const Label & GetAtomTypes() const
Access the atom types.
Definition: Snapshot.h:600
double wkj
Auxiliary variable w_kj.
Definition: COPSSImage.h:131
double xlo
Lower value for x.
Definition: COPSSImage.h:55
void force_pol(Snapshot *, size_t, size_t, size_t)
Image method kernel function, 1st order currently.
Definition: COPSSImage.cpp:134
double aa
Auxiliary variable.
Definition: COPSSImage.h:170
double vkj
Auxiliary variable v_kj.
Definition: COPSSImage.h:128
double Rxkj
Auxiliary variable R_kj for image kernel functions (x component).
Definition: COPSSImage.h:95
double e_
Gauss integration auxiliary parameter e.
Definition: COPSSImage.h:89
double aux3_integ
Auxiliary variable for integration.
Definition: COPSSImage.h:149
double xg0_[5]
Magic numbers for x.
Definition: COPSSImage.h:64
double Rij2
Auxiliary variable R_ij squared.
Definition: COPSSImage.h:119
double einner_
Dielectric constant of polarizable particles.
Definition: COPSSImage.h:43
double GetDielectric() const
Get dielectric constant.
Definition: Snapshot.h:169
double Getqqrd2e() const
Get qqrd2e value.
Definition: Snapshot.h:175
double auxv_y_delta
Auxiliary variable, delta (y component).
Definition: COPSSImage.h:158
double rkj
Auxiliary variable r_kj.
Definition: COPSSImage.h:107
double qqrd2e_
unit conversion constant.
Definition: COPSSImage.h:52
double aux3_delta
Auxiliary variable, delta.
Definition: COPSSImage.h:164
double aux1
Auxiliary variable.
Definition: COPSSImage.h:134
double Rzij
Auxiliary variable R_ij (z component).
Definition: COPSSImage.h:116
std::vector< double > atomTypeRadius_
List of radii for all atom types.
Definition: COPSSImage.h:70
double Ryij
Auxiliary variable R_ij (y component).
Definition: COPSSImage.h:113
double aux3Sqrt_integ
Auxiliary variable, square root of integration result.
Definition: COPSSImage.h:152
double auxv_y_integ
Auxiliary variable for integration (y component).
Definition: COPSSImage.h:143
int ngauss
Number of gaussian integrations.
Definition: COPSSImage.h:61
double auxv_z_delta
Auxiliary variable, delta (z component).
Definition: COPSSImage.h:161
double Rzkj
Auxiliary variable R_kj for image kernel functions (z component).
Definition: COPSSImage.h:101
double auxv_x_delta
Auxiliary variable, delta (x component).
Definition: COPSSImage.h:155
void PostSimulation(Snapshot *, const CVList &) override
Post-simulation hook.
Definition: COPSSImage.cpp:128
double auxv_z_integ
Auxiliary variable for integration (z component).
Definition: COPSSImage.h:146
double ginv_
Gauss integration auxiliary parameter g inverse.
Definition: COPSSImage.h:92
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:323
double Rkj2
Auxiliary variable R_kj squared.
Definition: COPSSImage.h:104
const std::vector< double > & GetCharges() const
Access the atom charges.
Definition: Snapshot.h:587
double aux2
Yet another auxiliary variable.
Definition: COPSSImage.h:137