mengine/src/ebufchrg.c

#define EXTERN extern

#include "pcwin.h"
#include "pcmod.h"

#include "energies.h"
#include "angles.h"
#include "derivs.h"
#include "hess.h"
#include "utility.h"
#include "field.h"
#include "attached.h"
#include "cutoffs.h"
#include "pot.h"
#include "solv.h"
                
EXTERN struct t_minim_values {
        int iprint, ndc, nconst;
        float dielc;
        } minim_values;
        
        
void ebufcharge()
{
   int i, j, k,temp;
   double xi, yi, zi, xk,yk,zk;
   double xr, yr, zr;
   double e, rik, cc, rik2;
   double rdielc;
   double rdn;
   double cutoff;

   cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
   
   energies.eu = 0.0;
   cc = 332.0716;
   rdielc = 1.0/units.dielec;

   rdn = 1.0;
   if (field.type == MMX || field.type == MM2 )
     rdn = .915;
   else if (field.type == MM3)
     rdn = .923;

   if (minim_values.iprint)
   {
       fprintf(pcmlogfile,"\nCharge-Charge Interactions\n");
       fprintf(pcmlogfile,"      At1      At2      q1        q2         Rik          Eqq\n");
   }

   for (i=1; i < natom; i++)
   {
      if (atom[i].charge != 0.0 )
      {
         
         for(j=i+1; j <= natom; j++)
         {
           if (atom[i].use || atom[j].use)
           {
             if ( skip[i][j] != i)
             {
               if (atom[j].charge != 0.0)
               { 
                  if (atom[i].type == 5 || atom[i].type == 36)
                  {
                   // compute reduced distance
                    xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
                    yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
                    zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
                  } else
                  {
                    xi = atom[i].x;
                    yi = atom[i].y;
                    zi = atom[i].z;
                  }
                  if (atom[j].type == 5 || atom[j].type == 36)
                  {
                   // compute reduced distance
                     xk = rdn*(atom[j].x-atom[atom[j].iat[0]].x) + atom[atom[j].iat[0]].x;
                     yk = rdn*(atom[j].y-atom[atom[j].iat[0]].y) + atom[atom[j].iat[0]].y;
                     zk = rdn*(atom[j].z-atom[atom[j].iat[0]].z) + atom[atom[j].iat[0]].z;
                  } else
                  {
                     xk = atom[j].x;
                     yk = atom[j].y;
                     zk = atom[j].z;
                  }
                  xr = xi - xk;
                  yr = yi - yk;
                  zr = zi - zk;
                  rik2 =  xr*xr + yr*yr + zr*zr;
                  if (rik2 < cutoff)
                  {
                     rik = sqrt( xr*xr + yr*yr + zr*zr);
                     e = cc*atom[i].charge*atom[j].charge/(units.dielec*(rik+0.05));
                     if (skip[i][j] == -i)
                       e *= units.chgscale;  
                     energies.eu += e;
                     atom[i].energy += e;
                     atom[j].energy += e;
                     if (minim_values.iprint)
                      fprintf(pcmlogfile,"QQ: %2s(%-3d)- %2s(%-3d)   %-8.3f  %-8.3f   %-8.3f = %8.4f\n",
                       atom[i].name, i, atom[j].name, j,atom[i].charge, atom[j].charge,rik,e);
                  }
                }
             }
           }
         }
      }
   }
}

void ebufcharge1()
{
   int i,j, k;
   double xi, yi, zi;
   double xr, yr, zr;
   double xk,yk,zk,rdn;
   double e, rik, cc, rik2, st;
   double rdielc;
   double cutoff;

   cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
   
   energies.eu = 0.0;
   cc = 332.0716;
   rdielc = 1.0/units.dielec;
   rdn = 1.0;
   if (field.type == MMX || field.type == MM2 )
     rdn = .915;
   else if (field.type == MM3)
     rdn = .923;

   for (i=0; i <= natom; i++)
   {
       deriv.deqq[i][0] = 0.0;
       deriv.deqq[i][1] = 0.0;
       deriv.deqq[i][2] = 0.0;
   }
      

   for (i=1; i < natom; i++)
   {
       if (atom[i].charge != 0.0 )
       {
           for(j=i+1; j <= natom; j++)
           {
             if (atom[i].use || atom[j].use)
             {
               if (atom[j].charge != 0.0 &&  skip[i][j] != i  )
               {
                if (atom[i].type == 5 || atom[i].type == 36)
                {
                   // compute reduced distance
                   xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
                   yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
                   zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
               } else
               {
                  xi = atom[i].x;
                  yi = atom[i].y;
                  zi = atom[i].z;
               }
                if (atom[j].type == 5 || atom[j].type == 36)
                {
                   // compute reduced distance
                   xk = rdn*(atom[j].x-atom[atom[j].iat[0]].x) + atom[atom[j].iat[0]].x;
                   yk = rdn*(atom[j].y-atom[atom[j].iat[0]].y) + atom[atom[j].iat[0]].y;
                   zk = rdn*(atom[j].z-atom[atom[j].iat[0]].z) + atom[atom[j].iat[0]].z;
               } else
               {
                  xk = atom[j].x;
                  yk = atom[j].y;
                  zk = atom[j].z;
               }
                xr = xi - xk;
                yr = yi - yk;
                zr = zi - zk;
                rik2 = xr*xr +yr*yr + zr*zr;
                if (rik2 < cutoff)
                {
                   rik = sqrt( rik2);
                   st = -cc*rdielc*atom[i].charge*atom[j].charge/(rik*(rik+0.05)*(rik+0.05));
                   if (skip[i][j] == -i)
                     st *= units.chgscale;
                   e = cc*rdielc*atom[i].charge*atom[j].charge/(rik+0.05);
                   if (skip[i][j] == -i)
                     e *= units.chgscale;
                  
                   deriv.deqq[i][0] += st*xr;
                   deriv.deqq[i][1] += st*yr;
                   deriv.deqq[i][2] += st*zr;
         
                   deriv.deqq[j][0] -= st*xr;
                   deriv.deqq[j][1] -= st*yr;
                   deriv.deqq[j][2] -= st*zr;
                            
                   energies.eu += e;
                }
               }
             }
           }
       }
   }
}

void ebufcharge2(int i)
{
    int j,k, jj,temp;
    double fik,de,d2e,d2edx,d2edy,d2edz;
    double xi,yi,zi,xr,yr,zr,term[3][3];
    double xk,yk,zk,rdn;
    double r,r2;
    double cc, rdielc;
    double cutoff;
    
    if (atom[i].charge == 0.0)
       return;

    cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
    
    cc = 332.0716;
    rdielc = 1.0/units.dielec;
    rdn = 1.0;
    if (field.type == MMX || field.type == MM2 )
     rdn = .915;
    else if (field.type == MM3)
     rdn = .923;
    
    
    if (atom[i].type == 5 || atom[i].type == 36)
    {
          // compute reduced distance
          xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
          yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
          zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
    } else
    {
         xi = atom[i].x;
         yi = atom[i].y;
         zi = atom[i].z;
     }
    
    for (k=1; k <= natom; k++)
    {
        if (atom[k].charge != 0.0  &&  (skip[i][k] != i) )
        {
                if (atom[j].type == 5 || atom[j].type == 36)
                {
                   // compute reduced distance
                   xk = rdn*(atom[k].x-atom[atom[k].iat[0]].x) + atom[atom[k].iat[0]].x;
                   yk = rdn*(atom[k].y-atom[atom[k].iat[0]].y) + atom[atom[k].iat[0]].y;
                   zk = rdn*(atom[k].z-atom[atom[k].iat[0]].z) + atom[atom[k].iat[0]].z;
               } else
               {
                  xk = atom[k].x;
                  yk = atom[k].y;
                  zk = atom[k].z;
               }
               xr = xi - xk;
               yr = yi - yk;
               zr = zi - zk;
              r2 = xr*xr + yr*yr + zr*zr;
              if (r2 < cutoff)
             {
              r = sqrt(r2);
              fik = cc*rdielc*atom[i].charge*atom[k].charge;
              if (skip[i][k] == -i)
                 fik *= units.chgscale;
              de = -fik/((r+0.05)*(r+0.05));
              d2e = -2.0*de/r;
              de = de / r;
              d2e = (d2e-de) / r2;
              d2edx = d2e * xr;
              d2edy = d2e * yr;
              d2edz = d2e * zr;
              term[0][0]= d2edx*xr + de;
              term[1][0] = d2edx*yr;
              term[2][0] = d2edx*zr;
              term[0][1] = term[1][0];
              term[1][1] = d2edy*yr + de;
              term[2][1] = d2edy*zr;
              term[0][2] = term[2][0];
              term[1][2] = term[2][1];
              term[2][2] = d2edz*zr + de;

              for (j=0; j < 3; j++)
              {
                  hess.hessx[i][j] += term[j][0];
                  hess.hessy[i][j] += term[j][1];
                  hess.hessz[i][j] += term[j][2];
                  hess.hessx[k][j] -= term[j][0];
                  hess.hessy[k][j] -= term[j][1];
                  hess.hessz[k][j] -= term[j][2];
              }
             }
        }
    }
}
Revision:	93
Committed:	Sat Jan 17 23:24:55 2009 UTC (13 years, 5 months ago) by wdelano
File size:	9245 byte(s)
Log Message:	branch created
Line	File contents
1	#define EXTERN extern
2
3	#include "pcwin.h"
4	#include "pcmod.h"
5
6	#include "energies.h"
7	#include "angles.h"
8	#include "derivs.h"
9	#include "hess.h"
10	#include "utility.h"
11	#include "field.h"
12	#include "attached.h"
13	#include "cutoffs.h"
14	#include "pot.h"
15	#include "solv.h"
16
17	EXTERN struct t_minim_values {
18	int iprint, ndc, nconst;
19	float dielc;
20	} minim_values;
21
22
23	void ebufcharge()
24	{
25	int i, j, k,temp;
26	double xi, yi, zi, xk,yk,zk;
27	double xr, yr, zr;
28	double e, rik, cc, rik2;
29	double rdielc;
30	double rdn;
31	double cutoff;
32
33	cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
34
35	energies.eu = 0.0;
36	cc = 332.0716;
37	rdielc = 1.0/units.dielec;
38
39	rdn = 1.0;
40	if (field.type == MMX \|\| field.type == MM2 )
41	rdn = .915;
42	else if (field.type == MM3)
43	rdn = .923;
44
45	if (minim_values.iprint)
46	{
47	fprintf(pcmlogfile,"\nCharge-Charge Interactions\n");
48	fprintf(pcmlogfile," At1 At2 q1 q2 Rik Eqq\n");
49	}
50
51	for (i=1; i < natom; i++)
52	{
53	if (atom[i].charge != 0.0 )
54	{
55
56	for(j=i+1; j <= natom; j++)
57	{
58	if (atom[i].use \|\| atom[j].use)
59	{
60	if ( skip[i][j] != i)
61	{
62	if (atom[j].charge != 0.0)
63	{
64	if (atom[i].type == 5 \|\| atom[i].type == 36)
65	{
66	// compute reduced distance
67	xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
68	yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
69	zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
70	} else
71	{
72	xi = atom[i].x;
73	yi = atom[i].y;
74	zi = atom[i].z;
75	}
76	if (atom[j].type == 5 \|\| atom[j].type == 36)
77	{
78	// compute reduced distance
79	xk = rdn*(atom[j].x-atom[atom[j].iat[0]].x) + atom[atom[j].iat[0]].x;
80	yk = rdn*(atom[j].y-atom[atom[j].iat[0]].y) + atom[atom[j].iat[0]].y;
81	zk = rdn*(atom[j].z-atom[atom[j].iat[0]].z) + atom[atom[j].iat[0]].z;
82	} else
83	{
84	xk = atom[j].x;
85	yk = atom[j].y;
86	zk = atom[j].z;
87	}
88	xr = xi - xk;
89	yr = yi - yk;
90	zr = zi - zk;
91	rik2 = xrxr + yryr + zr*zr;
92	if (rik2 < cutoff)
93	{
94	rik = sqrt( xrxr + yryr + zr*zr);
95	e = ccatom[i].chargeatom[j].charge/(units.dielec*(rik+0.05));
96	if (skip[i][j] == -i)
97	e *= units.chgscale;
98	energies.eu += e;
99	atom[i].energy += e;
100	atom[j].energy += e;
101	if (minim_values.iprint)
102	fprintf(pcmlogfile,"QQ: %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.3f = %8.4f\n",
103	atom[i].name, i, atom[j].name, j,atom[i].charge, atom[j].charge,rik,e);
104	}
105	}
106	}
107	}
108	}
109	}
110	}
111	}
112
113	void ebufcharge1()
114	{
115	int i,j, k;
116	double xi, yi, zi;
117	double xr, yr, zr;
118	double xk,yk,zk,rdn;
119	double e, rik, cc, rik2, st;
120	double rdielc;
121	double cutoff;
122
123	cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
124
125	energies.eu = 0.0;
126	cc = 332.0716;
127	rdielc = 1.0/units.dielec;
128	rdn = 1.0;
129	if (field.type == MMX \|\| field.type == MM2 )
130	rdn = .915;
131	else if (field.type == MM3)
132	rdn = .923;
133
134	for (i=0; i <= natom; i++)
135	{
136	deriv.deqq[i][0] = 0.0;
137	deriv.deqq[i][1] = 0.0;
138	deriv.deqq[i][2] = 0.0;
139	}
140
141
142	for (i=1; i < natom; i++)
143	{
144	if (atom[i].charge != 0.0 )
145	{
146	for(j=i+1; j <= natom; j++)
147	{
148	if (atom[i].use \|\| atom[j].use)
149	{
150	if (atom[j].charge != 0.0 && skip[i][j] != i )
151	{
152	if (atom[i].type == 5 \|\| atom[i].type == 36)
153	{
154	// compute reduced distance
155	xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
156	yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
157	zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
158	} else
159	{
160	xi = atom[i].x;
161	yi = atom[i].y;
162	zi = atom[i].z;
163	}
164	if (atom[j].type == 5 \|\| atom[j].type == 36)
165	{
166	// compute reduced distance
167	xk = rdn*(atom[j].x-atom[atom[j].iat[0]].x) + atom[atom[j].iat[0]].x;
168	yk = rdn*(atom[j].y-atom[atom[j].iat[0]].y) + atom[atom[j].iat[0]].y;
169	zk = rdn*(atom[j].z-atom[atom[j].iat[0]].z) + atom[atom[j].iat[0]].z;
170	} else
171	{
172	xk = atom[j].x;
173	yk = atom[j].y;
174	zk = atom[j].z;
175	}
176	xr = xi - xk;
177	yr = yi - yk;
178	zr = zi - zk;
179	rik2 = xrxr +yryr + zr*zr;
180	if (rik2 < cutoff)
181	{
182	rik = sqrt( rik2);
183	st = -ccrdielcatom[i].chargeatom[j].charge/(rik(rik+0.05)*(rik+0.05));
184	if (skip[i][j] == -i)
185	st *= units.chgscale;
186	e = ccrdielcatom[i].charge*atom[j].charge/(rik+0.05);
187	if (skip[i][j] == -i)
188	e *= units.chgscale;
189
190	deriv.deqq[i][0] += st*xr;
191	deriv.deqq[i][1] += st*yr;
192	deriv.deqq[i][2] += st*zr;
193
194	deriv.deqq[j][0] -= st*xr;
195	deriv.deqq[j][1] -= st*yr;
196	deriv.deqq[j][2] -= st*zr;
197
198	energies.eu += e;
199	}
200	}
201	}
202	}
203	}
204	}
205	}
206
207	void ebufcharge2(int i)
208	{
209	int j,k, jj,temp;
210	double fik,de,d2e,d2edx,d2edy,d2edz;
211	double xi,yi,zi,xr,yr,zr,term[3][3];
212	double xk,yk,zk,rdn;
213	double r,r2;
214	double cc, rdielc;
215	double cutoff;
216
217	if (atom[i].charge == 0.0)
218	return;
219
220	cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
221
222	cc = 332.0716;
223	rdielc = 1.0/units.dielec;
224	rdn = 1.0;
225	if (field.type == MMX \|\| field.type == MM2 )
226	rdn = .915;
227	else if (field.type == MM3)
228	rdn = .923;
229
230
231	if (atom[i].type == 5 \|\| atom[i].type == 36)
232	{
233	// compute reduced distance
234	xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
235	yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
236	zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
237	} else
238	{
239	xi = atom[i].x;
240	yi = atom[i].y;
241	zi = atom[i].z;
242	}
243
244	for (k=1; k <= natom; k++)
245	{
246	if (atom[k].charge != 0.0 && (skip[i][k] != i) )
247	{
248	if (atom[j].type == 5 \|\| atom[j].type == 36)
249	{
250	// compute reduced distance
251	xk = rdn*(atom[k].x-atom[atom[k].iat[0]].x) + atom[atom[k].iat[0]].x;
252	yk = rdn*(atom[k].y-atom[atom[k].iat[0]].y) + atom[atom[k].iat[0]].y;
253	zk = rdn*(atom[k].z-atom[atom[k].iat[0]].z) + atom[atom[k].iat[0]].z;
254	} else
255	{
256	xk = atom[k].x;
257	yk = atom[k].y;
258	zk = atom[k].z;
259	}
260	xr = xi - xk;
261	yr = yi - yk;
262	zr = zi - zk;
263	r2 = xrxr + yryr + zr*zr;
264	if (r2 < cutoff)
265	{
266	r = sqrt(r2);
267	fik = ccrdielcatom[i].charge*atom[k].charge;
268	if (skip[i][k] == -i)
269	fik *= units.chgscale;
270	de = -fik/((r+0.05)*(r+0.05));
271	d2e = -2.0*de/r;
272	de = de / r;
273	d2e = (d2e-de) / r2;
274	d2edx = d2e * xr;
275	d2edy = d2e * yr;
276	d2edz = d2e * zr;
277	term[0][0]= d2edx*xr + de;
278	term[1][0] = d2edx*yr;
279	term[2][0] = d2edx*zr;
280	term[0][1] = term[1][0];
281	term[1][1] = d2edy*yr + de;
282	term[2][1] = d2edy*zr;
283	term[0][2] = term[2][0];
284	term[1][2] = term[2][1];
285	term[2][2] = d2edz*zr + de;
286
287	for (j=0; j < 3; j++)
288	{
289	hess.hessx[i][j] += term[j][0];
290	hess.hessy[i][j] += term[j][1];
291	hess.hessz[i][j] += term[j][2];
292	hess.hessx[k][j] -= term[j][0];
293	hess.hessy[k][j] -= term[j][1];
294	hess.hessz[k][j] -= term[j][2];
295	}
296	}
297	}
298	}
299	}