?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Improved version for builds from modified trees.
[gromacs/qmmm-gamess-us.git] / src / mdlib / genborn.c
blob87848c18979cf6b3f51e77f3f1d0ccbb966a6eb9
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2008, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
26 *
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
34 */
35
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <math.h>
41 #include <string.h>
42
43 #include "typedefs.h"
44 #include "smalloc.h"
45 #include "genborn.h"
46 #include "vec.h"
47 #include "grompp.h"
48 #include "pdbio.h"
49 #include "names.h"
50 #include "physics.h"
51 #include "partdec.h"
52 #include "domdec.h"
53 #include "network.h"
54 #include "gmx_fatal.h"
55 #include "mtop_util.h"
56 #include "pbc.h"
57 #include "nrnb.h"
58 #include "bondf.h"
59
60 #ifdef GMX_LIB_MPI
61 #include <mpi.h>
62 #endif
63 #ifdef GMX_THREADS
64 #include "tmpi.h"
65 #endif
66
67 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_SSE2) )
68 #ifdef GMX_DOUBLE
69 #include "genborn_sse2_double.h"
70 #else
71 #include "genborn_sse2_single.h"
72 #include "genborn_allvsall_sse2_single.h"
73 #endif /* GMX_DOUBLE */
74 #endif /* GMX_SSE */
75
76 #include "genborn_allvsall.h"
77
78 //#define DISABLE_SSE
79
80 typedef struct {
81 int shift;
82 int naj;
83 int *aj;
84 int aj_nalloc;
85 } gbtmpnbl_t;
86
87 typedef struct gbtmpnbls {
88 int nlist;
89 gbtmpnbl_t *list;
90 int list_nalloc;
91 } t_gbtmpnbls;
92
93 /* This function is exactly the same as the one in bondfree.c. The reason
94 * it is copied here is that the bonded gb-interactions are evaluated
95 * not in calc_bonds, but rather in calc_gb_forces
96 */
97 static int pbc_rvec_sub(const t_pbc *pbc,const rvec xi,const rvec xj,rvec dx)
98 {
99 if (pbc) {
100 return pbc_dx_aiuc(pbc,xi,xj,dx);
101 }
102 else {
103 rvec_sub(xi,xj,dx);
104 return CENTRAL;
105 }
106 }
107
108 int init_gb_nblist(int natoms, t_nblist *nl)
109 {
110 nl->maxnri = natoms*4;
111 nl->maxnrj = 0;
112 nl->maxlen = 0;
113 nl->nri = 0;
114 nl->nrj = 0;
115 nl->iinr = NULL;
116 nl->gid = NULL;
117 nl->shift = NULL;
118 nl->jindex = NULL;
119 nl->jjnr = NULL;
120 /*nl->nltype = nltype;*/
121
122 srenew(nl->iinr, nl->maxnri);
123 srenew(nl->gid, nl->maxnri);
124 srenew(nl->shift, nl->maxnri);
125 srenew(nl->jindex, nl->maxnri+1);
126
127 nl->jindex[0] = 0;
128
129 return 0;
130 }
131
132 int print_nblist(int natoms, t_nblist *nl)
133 {
134 int i,k,ai,aj,nj0,nj1;
135
136 printf("genborn.c: print_nblist, natoms=%d\n",nl->nri);
137
138 for(i=0;i<nl->nri;i++)
139 {
140 ai=nl->iinr[i];
141 nj0=nl->jindex[i];
142 nj1=nl->jindex[i+1];
143
144 for(k=nj0;k<nj1;k++)
145 {
146 aj=nl->jjnr[k];
147 printf("ai=%d, aj=%d\n",ai,aj);
148 }
149 }
150
151 return 0;
152 }
153
154 typedef union {
155 real numlog;
156 int exp;
157 } u_table;
158
159 void fill_log_table(const int n, real *table)
160 {
161 u_table log_table;
162 real logfactor;
163 int i;
164
165 int incr = 1 << (23-n);
166 int p=pow(2,n);
167
168 logfactor = 1.0/log(2.0);
169
170 log_table.exp = 0x3F800000;
171
172 for(i=0;i<p;++i)
173 {
174 /* log2(numlog)=log(numlog)/log(2.0) */
175 table[i]=log(log_table.numlog)*logfactor;
176 log_table.exp+=incr;
177 }
178 }
179
180
181 real table_log(real val, const real *table, const int n)
182 {
183 int *const exp_ptr = ((int*)&val);
184 int x = *exp_ptr;
185 const int log_2 = ((x>>23) & 255) - 127;
186 x &= 0x7FFFFF;
187 x = x >> (23-n);
188 val = table[x];
189 return ((val+log_2)*0.69314718);
190 }
191
192 void gb_pd_send(t_commrec *cr, real *send_data, int nr)
193 {
194 #ifdef GMX_MPI
195 int i,cur;
196 int *index,*sendc,*disp;
197
198 snew(sendc,cr->nnodes);
199 snew(disp,cr->nnodes);
200
201 index = pd_index(cr);
202 cur = cr->nodeid;
203
204 /* Setup count/index arrays */
205 for(i=0;i<cr->nnodes;i++)
206 {
207 sendc[i] = index[i+1]-index[i];
208 disp[i] = index[i];
209 }
210
211 /* Do communication */
212 MPI_Gatherv(send_data+index[cur],sendc[cur],GMX_MPI_REAL,send_data,sendc,
213 disp,GMX_MPI_REAL,0,cr->mpi_comm_mygroup);
214 MPI_Bcast(send_data,nr,GMX_MPI_REAL,0,cr->mpi_comm_mygroup);
215
216 #endif
217 }
218
219
220 int init_gb_plist(t_params *p_list)
221 {
222 p_list->nr = 0;
223 p_list->param = NULL;
224
225 return 0;
226 }
227
228
229
230 int init_gb_still(const t_commrec *cr, t_forcerec *fr,
231 const t_atomtypes *atype, t_idef *idef, t_atoms *atoms,
232 gmx_genborn_t *born,int natoms)
233 {
234
235 int i,j,i1,i2,k,m,nbond,nang,ia,ib,ic,id,nb,idx,idx2,at;
236 int iam,ibm;
237 int at0,at1;
238 real length,angle;
239 real r,ri,rj,ri2,ri3,rj2,r2,r3,r4,rk,ratio,term,h,doffset;
240 real p1,p2,p3,factor,cosine,rab,rbc;
241
242 real *vsol;
243 real *gp;
244
245 snew(vsol,natoms);
246 snew(gp,natoms);
247 snew(born->gpol_still_work,natoms+3);
248
249 if(PAR(cr))
250 {
251 if(PARTDECOMP(cr))
252 {
253 pd_at_range(cr,&at0,&at1);
254
255 for(i=0;i<natoms;i++)
256 {
257 vsol[i] = gp[i] = 0;
258 }
259 }
260 else
261 {
262 at0 = 0;
263 at1 = natoms;
264 }
265 }
266 else
267 {
268 at0 = 0;
269 at1 = natoms;
270 }
271
272 doffset = born->gb_doffset;
273
274 for(i=0;i<natoms;i++)
275 {
276 born->gpol_globalindex[i]=born->vsolv_globalindex[i]=
277 born->gb_radius_globalindex[i]=0;
278 }
279
280 /* Compute atomic solvation volumes for Still method */
281 for(i=0;i<natoms;i++)
282 {
283 ri=atype->gb_radius[atoms->atom[i].type];
284 born->gb_radius_globalindex[i] = ri;
285 r3=ri*ri*ri;
286 born->vsolv_globalindex[i]=(4*M_PI/3)*r3;
287 }
288
289 for(j=0;j<idef->il[F_GB12].nr;j+=3)
290 {
291 m=idef->il[F_GB12].iatoms[j];
292 ia=idef->il[F_GB12].iatoms[j+1];
293 ib=idef->il[F_GB12].iatoms[j+2];
294
295 r=1.01*idef->iparams[m].gb.st;
296
297 ri = atype->gb_radius[atoms->atom[ia].type];
298 rj = atype->gb_radius[atoms->atom[ib].type];
299
300 ri2 = ri*ri;
301 ri3 = ri2*ri;
302 rj2 = rj*rj;
303
304 ratio = (rj2-ri2-r*r)/(2*ri*r);
305 h = ri*(1+ratio);
306 term = (M_PI/3.0)*h*h*(3.0*ri-h);
307
308 if(PARTDECOMP(cr))
309 {
310 vsol[ia]+=term;
311 }
312 else
313 {
314 born->vsolv_globalindex[ia] -= term;
315 }
316
317 ratio = (ri2-rj2-r*r)/(2*rj*r);
318 h = rj*(1+ratio);
319 term = (M_PI/3.0)*h*h*(3.0*rj-h);
320
321 if(PARTDECOMP(cr))
322 {
323 vsol[ib]+=term;
324 }
325 else
326 {
327 born->vsolv_globalindex[ib] -= term;
328 }
329 }
330
331 if(PARTDECOMP(cr))
332 {
333 gmx_sum(natoms,vsol,cr);
334
335 for(i=0;i<natoms;i++)
336 {
337 born->vsolv_globalindex[i]=born->vsolv_globalindex[i]-vsol[i];
338 }
339 }
340
341 /* Get the self-, 1-2 and 1-3 polarization energies for analytical Still
342 method */
343 /* Self */
344 for(j=0;j<natoms;j++)
345 {
346 if(born->use_globalindex[j]==1)
347 {
348 born->gpol_globalindex[j]=-0.5*ONE_4PI_EPS0/
349 (atype->gb_radius[atoms->atom[j].type]-doffset+STILL_P1);
350 }
351 }
352
353 /* 1-2 */
354 for(j=0;j<idef->il[F_GB12].nr;j+=3)
355 {
356 m=idef->il[F_GB12].iatoms[j];
357 ia=idef->il[F_GB12].iatoms[j+1];
358 ib=idef->il[F_GB12].iatoms[j+2];
359
360 r=idef->iparams[m].gb.st;
361
362 r4=r*r*r*r;
363
364 if(PARTDECOMP(cr))
365 {
366 gp[ia]+=STILL_P2*born->vsolv_globalindex[ib]/r4;
367 gp[ib]+=STILL_P2*born->vsolv_globalindex[ia]/r4;
368 }
369 else
370 {
371 born->gpol_globalindex[ia]=born->gpol_globalindex[ia]+
372 STILL_P2*born->vsolv_globalindex[ib]/r4;
373 born->gpol_globalindex[ib]=born->gpol_globalindex[ib]+
374 STILL_P2*born->vsolv_globalindex[ia]/r4;
375 }
376 }
377
378 /* 1-3 */
379 for(j=0;j<idef->il[F_GB13].nr;j+=3)
380 {
381 m=idef->il[F_GB13].iatoms[j];
382 ia=idef->il[F_GB13].iatoms[j+1];
383 ib=idef->il[F_GB13].iatoms[j+2];
384
385 r=idef->iparams[m].gb.st;
386 r4=r*r*r*r;
387
388 if(PARTDECOMP(cr))
389 {
390 gp[ia]+=STILL_P3*born->vsolv[ib]/r4;
391 gp[ib]+=STILL_P3*born->vsolv[ia]/r4;
392 }
393 else
394 {
395 born->gpol_globalindex[ia]=born->gpol_globalindex[ia]+
396 STILL_P3*born->vsolv_globalindex[ib]/r4;
397 born->gpol_globalindex[ib]=born->gpol_globalindex[ib]+
398 STILL_P3*born->vsolv_globalindex[ia]/r4;
399 }
400 }
401
402 if(PARTDECOMP(cr))
403 {
404 gmx_sum(natoms,gp,cr);
405
406 for(i=0;i<natoms;i++)
407 {
408 born->gpol_globalindex[i]=born->gpol_globalindex[i]+gp[i];
409 }
410 }
411
412 sfree(vsol);
413 sfree(gp);
414
415 return 0;
416 }
417
418
419
420 #define LOG_TABLE_ACCURACY 15 /* Accuracy of the table logarithm */
421
422
423 /* Initialize all GB datastructs and compute polarization energies */
424 int init_gb(gmx_genborn_t **p_born,
425 const t_commrec *cr, t_forcerec *fr, const t_inputrec *ir,
426 const gmx_mtop_t *mtop, real rgbradii, int gb_algorithm)
427 {
428 int i,j,m,ai,aj,jj,natoms,nalloc;
429 real rai,sk,p,doffset;
430
431 t_atoms atoms;
432 gmx_genborn_t *born;
433 gmx_localtop_t *localtop;
434
435 natoms = mtop->natoms;
436
437 atoms = gmx_mtop_global_atoms(mtop);
438 localtop = gmx_mtop_generate_local_top(mtop,ir);
439
440 snew(born,1);
441 *p_born = born;
442
443 born->nr = fr->natoms_force;
444 born->nr = natoms;
445
446 snew(born->drobc, natoms);
447 snew(born->bRad, natoms);
448
449 /* Allocate memory for the global data arrays */
450 snew(born->param_globalindex, natoms+3);
451 snew(born->gpol_globalindex, natoms+3);
452 snew(born->vsolv_globalindex, natoms+3);
453 snew(born->gb_radius_globalindex, natoms+3);
454 snew(born->use_globalindex, natoms+3);
455
456 snew(fr->invsqrta, natoms);
457 snew(fr->dvda, natoms);
458
459 fr->dadx = NULL;
460 fr->dadx_rawptr = NULL;
461 fr->nalloc_dadx = 0;
462 born->gpol_still_work = NULL;
463 born->gpol_hct_work = NULL;
464
465 /* snew(born->asurf,natoms); */
466 /* snew(born->dasurf,natoms); */
467
468 /* Initialize the gb neighbourlist */
469 init_gb_nblist(natoms,&(fr->gblist));
470
471 /* Do the Vsites exclusions (if any) */
472 for(i=0;i<natoms;i++)
473 {
474 jj = atoms.atom[i].type;
475 if (mtop->atomtypes.gb_radius[atoms.atom[i].type] > 0)
476 {
477 born->use_globalindex[i] = 1;
478 }
479 else
480 {
481 born->use_globalindex[i] = 0;
482 }
483
484 /* If we have a Vsite, put vs_globalindex[i]=0 */
485 if (C6 (fr->nbfp,fr->ntype,jj,jj) == 0 &&
486 C12(fr->nbfp,fr->ntype,jj,jj) == 0 &&
487 atoms.atom[i].q == 0)
488 {
489 born->use_globalindex[i]=0;
490 }
491 }
492
493 /* Copy algorithm parameters from inputrecord to local structure */
494 born->obc_alpha = ir->gb_obc_alpha;
495 born->obc_beta = ir->gb_obc_beta;
496 born->obc_gamma = ir->gb_obc_gamma;
497 born->gb_doffset = ir->gb_dielectric_offset;
498 born->gb_epsilon_solvent = ir->gb_epsilon_solvent;
499
500 doffset = born->gb_doffset;
501
502 /* If Still model, initialise the polarisation energies */
503 if(gb_algorithm==egbSTILL)
504 {
505 init_gb_still(cr, fr,&(mtop->atomtypes), &(localtop->idef), &atoms,
506 born, natoms);
507 }
508
509
510 /* If HCT/OBC, precalculate the sk*atype->S_hct factors */
511 else if(gb_algorithm==egbHCT || gb_algorithm==egbOBC)
512 {
513
514 snew(born->gpol_hct_work, natoms+3);
515
516 for(i=0;i<natoms;i++)
517 {
518 if(born->use_globalindex[i]==1)
519 {
520 rai = mtop->atomtypes.gb_radius[atoms.atom[i].type]-doffset;
521 sk = rai * mtop->atomtypes.S_hct[atoms.atom[i].type];
522 born->param_globalindex[i] = sk;
523 born->gb_radius_globalindex[i] = rai;
524 }
525 else
526 {
527 born->param_globalindex[i] = 0;
528 born->gb_radius_globalindex[i] = 0;
529 }
530 }
531 }
532
533 /* Init the logarithm table */
534 p=pow(2,LOG_TABLE_ACCURACY);
535 snew(born->log_table, p);
536
537 fill_log_table(LOG_TABLE_ACCURACY, born->log_table);
538
539 /* Allocate memory for work arrays for temporary use */
540 snew(born->work,natoms+4);
541 snew(born->count,natoms);
542 snew(born->nblist_work,natoms);
543
544 /* Domain decomposition specific stuff */
545 if(DOMAINDECOMP(cr))
546 {
547 snew(born->dd_work,natoms);
548 born->nlocal = cr->dd->nat_tot; /* cr->dd->nat_tot will be zero here */
549 }
550
551 return 0;
552 }
553
554
555
556 static int
557 calc_gb_rad_still(t_commrec *cr, t_forcerec *fr,int natoms, gmx_localtop_t *top,
558 const t_atomtypes *atype, rvec x[], t_nblist *nl,
559 gmx_genborn_t *born,t_mdatoms *md)
560 {
561 int i,k,n,nj0,nj1,ai,aj,type;
562 int shift;
563 real shX,shY,shZ;
564 real gpi,dr,dr2,dr4,idr4,rvdw,ratio,ccf,theta,term,rai,raj;
565 real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
566 real rinv,idr2,idr6,vaj,dccf,cosq,sinq,prod,gpi2;
567 real factor;
568 real vai, prod_ai, icf4,icf6;
569
570 factor = 0.5*ONE_4PI_EPS0;
571 n = 0;
572
573 for(i=0;i<born->nr;i++)
574 {
575 born->gpol_still_work[i]=0;
576 }
577
578 for(i=0;i<nl->nri;i++ )
579 {
580 ai = nl->iinr[i];
581
582 nj0 = nl->jindex[i];
583 nj1 = nl->jindex[i+1];
584
585 /* Load shifts for this list */
586 shift = nl->shift[i];
587 shX = fr->shift_vec[shift][0];
588 shY = fr->shift_vec[shift][1];
589 shZ = fr->shift_vec[shift][2];
590
591 gpi = 0;
592
593 rai = top->atomtypes.gb_radius[md->typeA[ai]];
594 vai = born->vsolv[ai];
595 prod_ai = STILL_P4*vai;
596
597 /* Load atom i coordinates, add shift vectors */
598 ix1 = shX + x[ai][0];
599 iy1 = shY + x[ai][1];
600 iz1 = shZ + x[ai][2];
601
602 for(k=nj0;k<nj1;k++)
603 {
604 aj = nl->jjnr[k];
605 jx1 = x[aj][0];
606 jy1 = x[aj][1];
607 jz1 = x[aj][2];
608
609 dx11 = ix1-jx1;
610 dy11 = iy1-jy1;
611 dz11 = iz1-jz1;
612
613 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
614 rinv = gmx_invsqrt(dr2);
615 idr2 = rinv*rinv;
616 idr4 = idr2*idr2;
617 idr6 = idr4*idr2;
618
619 raj = top->atomtypes.gb_radius[md->typeA[aj]];
620
621 rvdw = rai + raj;
622
623 ratio = dr2 / (rvdw * rvdw);
624 vaj = born->vsolv[aj];
625
626 if(ratio>STILL_P5INV)
627 {
628 ccf=1.0;
629 dccf=0.0;
630 }
631 else
632 {
633 theta = ratio*STILL_PIP5;
634 cosq = cos(theta);
635 term = 0.5*(1.0-cosq);
636 ccf = term*term;
637 sinq = 1.0 - cosq*cosq;
638 dccf = 2.0*term*sinq*gmx_invsqrt(sinq)*theta;
639 }
640
641 prod = STILL_P4*vaj;
642 icf4 = ccf*idr4;
643 icf6 = (4*ccf-dccf)*idr6;
644
645 born->gpol_still_work[aj] += prod_ai*icf4;
646 gpi = gpi+prod*icf4;
647
648 /* Save ai->aj and aj->ai chain rule terms */
649 fr->dadx[n++] = prod*icf6;
650 fr->dadx[n++] = prod_ai*icf6;
651 }
652 born->gpol_still_work[ai]+=gpi;
653 }
654
655 /* Parallel summations */
656 if(PARTDECOMP(cr))
657 {
658 gmx_sum(natoms, born->gpol_still_work, cr);
659 }
660 else if(DOMAINDECOMP(cr))
661 {
662 dd_atom_sum_real(cr->dd, born->gpol_still_work);
663 }
664
665 /* Calculate the radii */
666 for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
667 {
668 if(born->use[i] != 0)
669 {
670
671 gpi = born->gpol[i]+born->gpol_still_work[i];
672 gpi2 = gpi * gpi;
673 born->bRad[i] = factor*gmx_invsqrt(gpi2);
674 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
675 }
676 }
677
678 /* Extra communication required for DD */
679 if(DOMAINDECOMP(cr))
680 {
681 dd_atom_spread_real(cr->dd, born->bRad);
682 dd_atom_spread_real(cr->dd, fr->invsqrta);
683 }
684
685 return 0;
686
687 }
688
689
690 static int
691 calc_gb_rad_hct(t_commrec *cr,t_forcerec *fr,int natoms, gmx_localtop_t *top,
692 const t_atomtypes *atype, rvec x[], t_nblist *nl,
693 gmx_genborn_t *born,t_mdatoms *md)
694 {
695 int i,k,n,ai,aj,nj0,nj1,at0,at1;
696 int shift;
697 real shX,shY,shZ;
698 real rai,raj,gpi,dr2,dr,sk,sk_ai,sk2,sk2_ai,lij,uij,diff2,tmp,sum_ai;
699 real rad,min_rad,rinv,rai_inv;
700 real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
701 real lij2, uij2, lij3, uij3, t1,t2,t3;
702 real lij_inv,dlij,duij,sk2_rinv,prod,log_term;
703 real doffset,raj_inv,dadx_val;
704 real *gb_radius;
705
706 doffset = born->gb_doffset;
707 gb_radius = born->gb_radius;
708
709 for(i=0;i<born->nr;i++)
710 {
711 born->gpol_hct_work[i] = 0;
712 }
713
714 /* Keep the compiler happy */
715 n = 0;
716 prod = 0;
717
718 for(i=0;i<nl->nri;i++)
719 {
720 ai = nl->iinr[i];
721
722 nj0 = nl->jindex[ai];
723 nj1 = nl->jindex[ai+1];
724
725 /* Load shifts for this list */
726 shift = nl->shift[i];
727 shX = fr->shift_vec[shift][0];
728 shY = fr->shift_vec[shift][1];
729 shZ = fr->shift_vec[shift][2];
730
731 rai = gb_radius[ai];
732 rai_inv = 1.0/rai;
733
734 sk_ai = born->param[ai];
735 sk2_ai = sk_ai*sk_ai;
736
737 /* Load atom i coordinates, add shift vectors */
738 ix1 = shX + x[ai][0];
739 iy1 = shY + x[ai][1];
740 iz1 = shZ + x[ai][2];
741
742 sum_ai = 0;
743
744 for(k=nj0;k<nj1;k++)
745 {
746 aj = nl->jjnr[k];
747
748 jx1 = x[aj][0];
749 jy1 = x[aj][1];
750 jz1 = x[aj][2];
751
752 dx11 = ix1 - jx1;
753 dy11 = iy1 - jy1;
754 dz11 = iz1 - jz1;
755
756 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
757 rinv = gmx_invsqrt(dr2);
758 dr = rinv*dr2;
759
760 sk = born->param[aj];
761 raj = gb_radius[aj];
762
763 /* aj -> ai interaction */
764 if(rai < dr+sk)
765 {
766 lij = 1.0/(dr-sk);
767 dlij = 1.0;
768
769 if(rai>dr-sk)
770 {
771 lij = rai_inv;
772 dlij = 0.0;
773 }
774
775 lij2 = lij*lij;
776 lij3 = lij2*lij;
777
778 uij = 1.0/(dr+sk);
779 uij2 = uij*uij;
780 uij3 = uij2*uij;
781
782 diff2 = uij2-lij2;
783
784 lij_inv = gmx_invsqrt(lij2);
785 sk2 = sk*sk;
786 sk2_rinv = sk2*rinv;
787 prod = 0.25*sk2_rinv;
788
789 /* log_term = table_log(uij*lij_inv,born->log_table,
790 LOG_TABLE_ACCURACY); */
791 log_term = log(uij*lij_inv);
792
793 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
794 prod*(-diff2);
795
796 if(rai<sk-dr)
797 {
798 tmp = tmp + 2.0 * (rai_inv-lij);
799 }
800
801 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
802 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
803 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
804
805 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
806 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */
807 /* rb2 is moved to chainrule */
808
809 sum_ai += 0.5*tmp;
810 }
811 else
812 {
813 dadx_val = 0.0;
814 }
815 fr->dadx[n++] = dadx_val;
816
817
818 /* ai -> aj interaction */
819 if(raj < dr + sk_ai)
820 {
821 lij = 1.0/(dr-sk_ai);
822 dlij = 1.0;
823 raj_inv = 1.0/raj;
824
825 if(raj>dr-sk_ai)
826 {
827 lij = raj_inv;
828 dlij = 0.0;
829 }
830
831 lij2 = lij * lij;
832 lij3 = lij2 * lij;
833
834 uij = 1.0/(dr+sk_ai);
835 uij2 = uij * uij;
836 uij3 = uij2 * uij;
837
838 diff2 = uij2-lij2;
839
840 lij_inv = gmx_invsqrt(lij2);
841 sk2 = sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
842 sk2_rinv = sk2*rinv;
843 prod = 0.25 * sk2_rinv;
844
845 /* log_term = table_log(uij*lij_inv,born->log_table,
846 LOG_TABLE_ACCURACY); */
847 log_term = log(uij*lij_inv);
848
849 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
850 prod*(-diff2);
851
852 if(raj<sk_ai-dr)
853 {
854 tmp = tmp + 2.0 * (raj_inv-lij);
855 }
856
857 /* duij = 1.0 */
858 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
859 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
860 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
861
862 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
863 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */ /* rb2 is moved to chainrule */
864
865 born->gpol_hct_work[aj] += 0.5*tmp;
866 }
867 else
868 {
869 dadx_val = 0.0;
870 }
871 fr->dadx[n++] = dadx_val;
872 }
873
874 born->gpol_hct_work[ai] += sum_ai;
875 }
876
877 /* Parallel summations */
878 if(PARTDECOMP(cr))
879 {
880 gmx_sum(natoms, born->gpol_hct_work, cr);
881 }
882 else if(DOMAINDECOMP(cr))
883 {
884 dd_atom_sum_real(cr->dd, born->gpol_hct_work);
885 }
886
887 for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
888 {
889 if(born->use[i] != 0)
890 {
891 rai = top->atomtypes.gb_radius[md->typeA[i]]-doffset;
892 sum_ai = 1.0/rai - born->gpol_hct_work[i];
893 min_rad = rai + doffset;
894 rad = 1.0/sum_ai;
895
896 born->bRad[i] = rad > min_rad ? rad : min_rad;
897 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
898 }
899 }
900
901 /* Extra communication required for DD */
902 if(DOMAINDECOMP(cr))
903 {
904 dd_atom_spread_real(cr->dd, born->bRad);
905 dd_atom_spread_real(cr->dd, fr->invsqrta);
906 }
907
908
909 return 0;
910 }
911
912 static int
913 calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, int natoms, gmx_localtop_t *top,
914 const t_atomtypes *atype, rvec x[], t_nblist *nl, gmx_genborn_t *born,t_mdatoms *md)
915 {
916 int i,k,ai,aj,nj0,nj1,n,at0,at1;
917 int shift;
918 real shX,shY,shZ;
919 real rai,raj,gpi,dr2,dr,sk,sk2,lij,uij,diff2,tmp,sum_ai;
920 real rad, min_rad,sum_ai2,sum_ai3,tsum,tchain,rinv,rai_inv,lij_inv,rai_inv2;
921 real log_term,prod,sk2_rinv,sk_ai,sk2_ai;
922 real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
923 real lij2,uij2,lij3,uij3,dlij,duij,t1,t2,t3;
924 real doffset,raj_inv,dadx_val;
925 real *gb_radius;
926
927 /* Keep the compiler happy */
928 n = 0;
929 prod = 0;
930 raj = 0;
931
932 doffset = born->gb_doffset;
933 gb_radius = born->gb_radius;
934
935 for(i=0;i<born->nr;i++)
936 {
937 born->gpol_hct_work[i] = 0;
938 }
939
940 for(i=0;i<nl->nri;i++)
941 {
942 ai = nl->iinr[i];
943
944 nj0 = nl->jindex[i];
945 nj1 = nl->jindex[i+1];
946
947 /* Load shifts for this list */
948 shift = nl->shift[i];
949 shX = fr->shift_vec[shift][0];
950 shY = fr->shift_vec[shift][1];
951 shZ = fr->shift_vec[shift][2];
952
953 rai = gb_radius[ai];
954 rai_inv = 1.0/rai;
955
956 sk_ai = born->param[ai];
957 sk2_ai = sk_ai*sk_ai;
958
959 /* Load atom i coordinates, add shift vectors */
960 ix1 = shX + x[ai][0];
961 iy1 = shY + x[ai][1];
962 iz1 = shZ + x[ai][2];
963
964 sum_ai = 0;
965
966 for(k=nj0;k<nj1;k++)
967 {
968 aj = nl->jjnr[k];
969
970 jx1 = x[aj][0];
971 jy1 = x[aj][1];
972 jz1 = x[aj][2];
973
974 dx11 = ix1 - jx1;
975 dy11 = iy1 - jy1;
976 dz11 = iz1 - jz1;
977
978 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
979 rinv = gmx_invsqrt(dr2);
980 dr = dr2*rinv;
981
982 /* sk is precalculated in init_gb() */
983 sk = born->param[aj];
984 raj = gb_radius[aj];
985
986 /* aj -> ai interaction */
987 if(rai < dr+sk)
988 {
989 lij = 1.0/(dr-sk);
990 dlij = 1.0;
991
992 if(rai>dr-sk)
993 {
994 lij = rai_inv;
995 dlij = 0.0;
996 }
997
998 uij = 1.0/(dr+sk);
999 lij2 = lij * lij;
1000 lij3 = lij2 * lij;
1001 uij2 = uij * uij;
1002 uij3 = uij2 * uij;
1003
1004 diff2 = uij2-lij2;
1005
1006 lij_inv = gmx_invsqrt(lij2);
1007 sk2 = sk*sk;
1008 sk2_rinv = sk2*rinv;
1009 prod = 0.25*sk2_rinv;
1010
1011 log_term = log(uij*lij_inv);
1012
1013 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
1014 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
1015
1016 if(rai < sk-dr)
1017 {
1018 tmp = tmp + 2.0 * (rai_inv-lij);
1019 }
1020
1021 /* duij = 1.0; */
1022 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
1023 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
1024 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
1025
1026 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
1027
1028 sum_ai += 0.5*tmp;
1029 }
1030 else
1031 {
1032 dadx_val = 0.0;
1033 }
1034 fr->dadx[n++] = dadx_val;
1035
1036 /* ai -> aj interaction */
1037 if(raj < dr + sk_ai)
1038 {
1039 lij = 1.0/(dr-sk_ai);
1040 dlij = 1.0;
1041 raj_inv = 1.0/raj;
1042
1043 if(raj>dr-sk_ai)
1044 {
1045 lij = raj_inv;
1046 dlij = 0.0;
1047 }
1048
1049 lij2 = lij * lij;
1050 lij3 = lij2 * lij;
1051
1052 uij = 1.0/(dr+sk_ai);
1053 uij2 = uij * uij;
1054 uij3 = uij2 * uij;
1055
1056 diff2 = uij2-lij2;
1057
1058 lij_inv = gmx_invsqrt(lij2);
1059 sk2 = sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
1060 sk2_rinv = sk2*rinv;
1061 prod = 0.25 * sk2_rinv;
1062
1063 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
1064 log_term = log(uij*lij_inv);
1065
1066 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
1067
1068 if(raj<sk_ai-dr)
1069 {
1070 tmp = tmp + 2.0 * (raj_inv-lij);
1071 }
1072
1073 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
1074 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
1075 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
1076
1077 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
1078
1079 born->gpol_hct_work[aj] += 0.5*tmp;
1080
1081 }
1082 else
1083 {
1084 dadx_val = 0.0;
1085 }
1086 fr->dadx[n++] = dadx_val;
1087
1088 }
1089 born->gpol_hct_work[ai] += sum_ai;
1090
1091 }
1092
1093 /* Parallel summations */
1094 if(PARTDECOMP(cr))
1095 {
1096 gmx_sum(natoms, born->gpol_hct_work, cr);
1097 }
1098 else if(DOMAINDECOMP(cr))
1099 {
1100 dd_atom_sum_real(cr->dd, born->gpol_hct_work);
1101 }
1102
1103 for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
1104 {
1105 if(born->use[i] != 0)
1106 {
1107 rai = top->atomtypes.gb_radius[md->typeA[i]];
1108 rai_inv2 = 1.0/rai;
1109 rai = rai-doffset;
1110 rai_inv = 1.0/rai;
1111 sum_ai = rai * born->gpol_hct_work[i];
1112 sum_ai2 = sum_ai * sum_ai;
1113 sum_ai3 = sum_ai2 * sum_ai;
1114
1115 tsum = tanh(born->obc_alpha*sum_ai-born->obc_beta*sum_ai2+born->obc_gamma*sum_ai3);
1116 born->bRad[i] = rai_inv - tsum*rai_inv2;
1117 born->bRad[i] = 1.0 / born->bRad[i];
1118
1119 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
1120
1121 tchain = rai * (born->obc_alpha-2*born->obc_beta*sum_ai+3*born->obc_gamma*sum_ai2);
1122 born->drobc[i] = (1.0-tsum*tsum)*tchain*rai_inv2;
1123 }
1124 }
1125
1126 /* Extra (local) communication required for DD */
1127 if(DOMAINDECOMP(cr))
1128 {
1129 dd_atom_spread_real(cr->dd, born->bRad);
1130 dd_atom_spread_real(cr->dd, fr->invsqrta);
1131 dd_atom_spread_real(cr->dd, born->drobc);
1132 }
1133
1134 return 0;
1135
1136 }
1137
1138
1139
1140 int calc_gb_rad(t_commrec *cr, t_forcerec *fr, t_inputrec *ir,gmx_localtop_t *top,
1141 const t_atomtypes *atype, rvec x[], t_nblist *nl, gmx_genborn_t *born,t_mdatoms *md,t_nrnb *nrnb)
1142 {
1143 real *p;
1144 int cnt;
1145 int ndadx;
1146
1147 if(fr->bAllvsAll && fr->dadx==NULL)
1148 {
1149 /* We might need up to 8 atoms of padding before and after,
1150 * and another 4 units to guarantee SSE alignment.
1151 */
1152 fr->nalloc_dadx = 2*(md->homenr+12)*(md->nr/2+1+12);
1153 snew(fr->dadx_rawptr,fr->nalloc_dadx);
1154 fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
1155 }
1156 else
1157 {
1158 /* In the SSE-enabled gb-loops, when writing to dadx, we
1159 * always write 2*4 elements at a time, even in the case with only
1160 * 1-3 j particles, where we only really need to write 2*(1-3)
1161 * elements. This is because we want dadx to be aligned to a 16-
1162 * byte boundary, and being able to use _mm_store/load_ps
1163 */
1164 ndadx = 2 * (nl->nrj + 3*nl->nri);
1165
1166 /* First, reallocate the dadx array, we need 3 extra for SSE */
1167 if (ndadx + 3 > fr->nalloc_dadx)
1168 {
1169 fr->nalloc_dadx = over_alloc_large(ndadx) + 3;
1170 srenew(fr->dadx_rawptr,fr->nalloc_dadx);
1171 fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
1172 }
1173 }
1174
1175 #ifndef DOUBLE
1176 if(fr->bAllvsAll)
1177 {
1178 cnt = md->homenr*(md->nr/2+1);
1179
1180 if(ir->gb_algorithm==egbSTILL)
1181 {
1182 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_SSE2) )
1183 genborn_allvsall_calc_still_radii_sse2_single(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1184 #else
1185 genborn_allvsall_calc_still_radii(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1186 #endif
1187 inc_nrnb(nrnb,eNR_BORN_AVA_RADII_STILL,cnt);
1188 }
1189 else if(ir->gb_algorithm==egbHCT || ir->gb_algorithm==egbOBC)
1190 {
1191 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_SSE2) )
1192 genborn_allvsall_calc_hct_obc_radii_sse2_single(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1193 #else
1194 genborn_allvsall_calc_hct_obc_radii(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1195 #endif
1196 inc_nrnb(nrnb,eNR_BORN_AVA_RADII_HCT_OBC,cnt);
1197 }
1198 else
1199 {
1200 gmx_fatal(FARGS,"Bad gb algorithm for all-vs-all interactions");
1201 }
1202 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,md->homenr);
1203
1204 return 0;
1205 }
1206 #endif
1207
1208 /* Switch for determining which algorithm to use for Born radii calculation */
1209 #ifdef GMX_DOUBLE
1210
1211 #if ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || defined(GMX_SSE2) )
1212 /* x86 or x86-64 with GCC inline assembly and/or SSE intrinsics */
1213 switch(ir->gb_algorithm)
1214 {
1215 case egbSTILL:
1216 calc_gb_rad_still_sse2_double(cr,fr,md->nr,top, atype, x[0], nl, born, md);
1217 break;
1218 case egbHCT:
1219 calc_gb_rad_hct_sse2_double(cr,fr,md->nr,top, atype, x[0], nl, born, md);
1220 break;
1221 case egbOBC:
1222 calc_gb_rad_obc_sse2_double(cr,fr,md->nr,top, atype, x[0], nl, born, md);
1223 break;
1224
1225 default:
1226 gmx_fatal(FARGS, "Unknown double precision sse-enabled algorithm for Born radii calculation: %d",ir->gb_algorithm);
1227 }
1228 #else
1229 switch(ir->gb_algorithm)
1230 {
1231 case egbSTILL:
1232 calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1233 break;
1234 case egbHCT:
1235 calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1236 break;
1237 case egbOBC:
1238 calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1239 break;
1240
1241 default:
1242 gmx_fatal(FARGS, "Unknown double precision algorithm for Born radii calculation: %d",ir->gb_algorithm);
1243 }
1244
1245 #endif
1246
1247 #else
1248
1249 #if (!defined DISABLE_SSE && ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_SSE2) ) )
1250 /* x86 or x86-64 with GCC inline assembly and/or SSE intrinsics */
1251 switch(ir->gb_algorithm)
1252 {
1253 case egbSTILL:
1254 calc_gb_rad_still_sse(cr,fr,born->nr,top, atype, x[0], nl, born, md);
1255 break;
1256 case egbHCT:
1257 case egbOBC:
1258 calc_gb_rad_hct_obc_sse(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1259 break;
1260
1261 default:
1262 gmx_fatal(FARGS, "Unknown sse-enabled algorithm for Born radii calculation: %d",ir->gb_algorithm);
1263 }
1264
1265 #else
1266 switch(ir->gb_algorithm)
1267 {
1268 case egbSTILL:
1269 calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1270 break;
1271 case egbHCT:
1272 calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1273 break;
1274 case egbOBC:
1275 calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1276 break;
1277
1278 default:
1279 gmx_fatal(FARGS, "Unknown algorithm for Born radii calculation: %d",ir->gb_algorithm);
1280 }
1281
1282 #endif /* Single precision sse */
1283
1284 #endif /* Double or single precision */
1285
1286 if(fr->bAllvsAll==FALSE)
1287 {
1288 switch(ir->gb_algorithm)
1289 {
1290 case egbSTILL:
1291 inc_nrnb(nrnb,eNR_BORN_RADII_STILL,nl->nrj);
1292 break;
1293 case egbHCT:
1294 case egbOBC:
1295 inc_nrnb(nrnb,eNR_BORN_RADII_HCT_OBC,nl->nrj);
1296 break;
1297
1298 default:
1299 break;
1300 }
1301 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,nl->nri);
1302 }
1303
1304 return 0;
1305 }
1306
1307
1308
1309 real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
1310 real *invsqrta, real *dvda, real *GBtab, t_idef *idef,
1311 real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph)
1312 {
1313 int i,j,n0,m,nnn,type,ai,aj;
1314 int ki;
1315
1316 real isai,isaj;
1317 real r,rsq11;
1318 real rinv11,iq;
1319 real isaprod,qq,gbscale,gbtabscale,Y,F,Geps,Heps2,Fp,VV,FF,rt,eps,eps2;
1320 real vgb,fgb,vcoul,fijC,dvdatmp,fscal,dvdaj;
1321 real vctot;
1322
1323 rvec dx;
1324 ivec dt;
1325
1326 t_iatom *forceatoms;
1327
1328 /* Scale the electrostatics by gb_epsilon_solvent */
1329 facel = facel * (1.0 - 1.0/gb_epsilon_solvent);
1330
1331 gbtabscale=*p_gbtabscale;
1332 vctot = 0.0;
1333
1334 for(j=F_GB12;j<=F_GB14;j++)
1335 {
1336 forceatoms = idef->il[j].iatoms;
1337
1338 for(i=0;i<idef->il[j].nr; )
1339 {
1340 /* To avoid reading in the interaction type, we just increment i to pass over
1341 * the types in the forceatoms array, this saves some memory accesses
1342 */
1343 i++;
1344 ai = forceatoms[i++];
1345 aj = forceatoms[i++];
1346
1347 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx);
1348 rsq11 = iprod(dx,dx);
1349
1350 isai = invsqrta[ai];
1351 iq = (-1)*facel*charge[ai];
1352
1353 rinv11 = gmx_invsqrt(rsq11);
1354 isaj = invsqrta[aj];
1355 isaprod = isai*isaj;
1356 qq = isaprod*iq*charge[aj];
1357 gbscale = isaprod*gbtabscale;
1358 r = rsq11*rinv11;
1359 rt = r*gbscale;
1360 n0 = rt;
1361 eps = rt-n0;
1362 eps2 = eps*eps;
1363 nnn = 4*n0;
1364 Y = GBtab[nnn];
1365 F = GBtab[nnn+1];
1366 Geps = eps*GBtab[nnn+2];
1367 Heps2 = eps2*GBtab[nnn+3];
1368 Fp = F+Geps+Heps2;
1369 VV = Y+eps*Fp;
1370 FF = Fp+Geps+2.0*Heps2;
1371 vgb = qq*VV;
1372 fijC = qq*FF*gbscale;
1373 dvdatmp = -(vgb+fijC*r)*0.5;
1374 dvda[aj] = dvda[aj] + dvdatmp*isaj*isaj;
1375 dvda[ai] = dvda[ai] + dvdatmp*isai*isai;
1376 vctot = vctot + vgb;
1377 fgb = -(fijC)*rinv11;
1378
1379 if (graph) {
1380 ivec_sub(SHIFT_IVEC(graph,ai),SHIFT_IVEC(graph,aj),dt);
1381 ki=IVEC2IS(dt);
1382 }
1383
1384 for (m=0; (m<DIM); m++) { /* 15 */
1385 fscal=fgb*dx[m];
1386 f[ai][m]+=fscal;
1387 f[aj][m]-=fscal;
1388 fshift[ki][m]+=fscal;
1389 fshift[CENTRAL][m]-=fscal;
1390 }
1391 }
1392 }
1393
1394 return vctot;
1395 }
1396
1397 real calc_gb_selfcorrections(t_commrec *cr, int natoms,
1398 real *charge, gmx_genborn_t *born, real *dvda, t_mdatoms *md, double facel)
1399 {
1400 int i,ai,at0,at1;
1401 real rai,e,derb,q,q2,fi,rai_inv,vtot;
1402
1403 if(PARTDECOMP(cr))
1404 {
1405 pd_at_range(cr,&at0,&at1);
1406 }
1407 else if(DOMAINDECOMP(cr))
1408 {
1409 at0=0;
1410 at1=cr->dd->nat_home;
1411 }
1412 else
1413 {
1414 at0=0;
1415 at1=natoms;
1416
1417 }
1418
1419 /* Scale the electrostatics by gb_epsilon_solvent */
1420 facel = facel * (1.0 - 1.0/born->gb_epsilon_solvent);
1421
1422 vtot=0.0;
1423
1424 /* Apply self corrections */
1425 for(i=at0;i<at1;i++)
1426 {
1427 ai = i;
1428
1429 if(born->use[ai]==1)
1430 {
1431 rai = born->bRad[ai];
1432 rai_inv = 1.0/rai;
1433 q = charge[ai];
1434 q2 = q*q;
1435 fi = facel*q2;
1436 e = fi*rai_inv;
1437 derb = 0.5*e*rai_inv*rai_inv;
1438 dvda[ai] += derb*rai;
1439 vtot -= 0.5*e;
1440 }
1441 }
1442
1443 return vtot;
1444
1445 }
1446
1447 real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr,int natoms,gmx_genborn_t *born, gmx_localtop_t *top,
1448 const t_atomtypes *atype, real *dvda,int gb_algorithm, t_mdatoms *md)
1449 {
1450 int ai,i,at0,at1;
1451 real e,es,rai,rbi,term,probe,tmp,factor;
1452 real rbi_inv,rbi_inv2;
1453
1454 /* To keep the compiler happy */
1455 factor=0;
1456
1457 if(PARTDECOMP(cr))
1458 {
1459 pd_at_range(cr,&at0,&at1);
1460 }
1461 else if(DOMAINDECOMP(cr))
1462 {
1463 at0 = 0;
1464 at1 = cr->dd->nat_home;
1465 }
1466 else
1467 {
1468 at0=0;
1469 at1=natoms;
1470 }
1471
1472 /* The surface area factor is 0.0049 for Still model, 0.0054 for HCT/OBC */
1473 if(gb_algorithm==egbSTILL)
1474 {
1475 factor=0.0049*100*CAL2JOULE;
1476 }
1477 else
1478 {
1479 factor=0.0054*100*CAL2JOULE;
1480 }
1481
1482 /* if(gb_algorithm==egbHCT || gb_algorithm==egbOBC) */
1483
1484 es = 0;
1485 probe = 0.14;
1486 term = M_PI*4;
1487
1488 for(i=at0;i<at1;i++)
1489 {
1490 ai = i;
1491
1492 if(born->use[ai]==1)
1493 {
1494 rai = top->atomtypes.gb_radius[md->typeA[ai]];
1495 rbi_inv = fr->invsqrta[ai];
1496 rbi_inv2 = rbi_inv * rbi_inv;
1497 tmp = (rai*rbi_inv2)*(rai*rbi_inv2);
1498 tmp = tmp*tmp*tmp;
1499 e = factor*term*(rai+probe)*(rai+probe)*tmp;
1500 dvda[ai] = dvda[ai] - 6*e*rbi_inv2;
1501 es = es + e;
1502 }
1503 }
1504
1505 return es;
1506 }
1507
1508
1509
1510 real calc_gb_chainrule(int natoms, t_nblist *nl, real *dadx, real *dvda, rvec x[], rvec t[], rvec fshift[],
1511 rvec shift_vec[], int gb_algorithm, gmx_genborn_t *born, t_mdatoms *md)
1512 {
1513 int i,k,n,ai,aj,nj0,nj1,n0,n1;
1514 int shift;
1515 real shX,shY,shZ;
1516 real fgb,fij,rb2,rbi,fix1,fiy1,fiz1;
1517 real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11,rsq11;
1518 real rinv11,tx,ty,tz,rbai,rbaj,fgb_ai;
1519 real *rb;
1520 volatile int idx;
1521
1522 n = 0;
1523 rb = born->work;
1524
1525
1526 n0 = md->start;
1527 n1 = md->start+md->homenr+1+natoms/2;
1528
1529 if(gb_algorithm==egbSTILL)
1530 {
1531 for(i=n0;i<n1;i++)
1532 {
1533 k = i % natoms;
1534 rbi = born->bRad[k];
1535 rb[k] = (2 * rbi * rbi * dvda[k])/ONE_4PI_EPS0;
1536 }
1537 }
1538 else if(gb_algorithm==egbHCT)
1539 {
1540 for(i=n0;i<n1;i++)
1541 {
1542 k = i % natoms;
1543 rbi = born->bRad[k];
1544 rb[k] = rbi * rbi * dvda[k];
1545 }
1546 }
1547 else if(gb_algorithm==egbOBC)
1548 {
1549 for(i=n0;i<n1;i++)
1550 {
1551 k = i % natoms;
1552 rbi = born->bRad[k];
1553 rb[k] = rbi * rbi * born->drobc[k] * dvda[k];
1554 }
1555 }
1556
1557 for(i=0;i<nl->nri;i++)
1558 {
1559 ai = nl->iinr[i];
1560
1561 nj0 = nl->jindex[ai];
1562 nj1 = nl->jindex[ai+1];
1563
1564 /* Load shifts for this list */
1565 shift = nl->shift[i];
1566 shX = shift_vec[shift][0];
1567 shY = shift_vec[shift][1];
1568 shZ = shift_vec[shift][2];
1569
1570 /* Load atom i coordinates, add shift vectors */
1571 ix1 = shX + x[ai][0];
1572 iy1 = shY + x[ai][1];
1573 iz1 = shZ + x[ai][2];
1574
1575 fix1 = 0;
1576 fiy1 = 0;
1577 fiz1 = 0;
1578
1579 rbai = rb[ai];
1580
1581 for(k=nj0;k<nj1;k++)
1582 {
1583 aj = nl->jjnr[k];
1584
1585 jx1 = x[aj][0];
1586 jy1 = x[aj][1];
1587 jz1 = x[aj][2];
1588
1589 dx11 = ix1 - jx1;
1590 dy11 = iy1 - jy1;
1591 dz11 = iz1 - jz1;
1592
1593 rbaj = rb[aj];
1594
1595 fgb = rbai*dadx[n++];
1596 fgb_ai = rbaj*dadx[n++];
1597
1598 /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1599 fgb = fgb + fgb_ai;
1600
1601 tx = fgb * dx11;
1602 ty = fgb * dy11;
1603 tz = fgb * dz11;
1604
1605 fix1 = fix1 + tx;
1606 fiy1 = fiy1 + ty;
1607 fiz1 = fiz1 + tz;
1608
1609 /* Update force on atom aj */
1610 t[aj][0] = t[aj][0] - tx;
1611 t[aj][1] = t[aj][1] - ty;
1612 t[aj][2] = t[aj][2] - tz;
1613 }
1614
1615 /* Update force and shift forces on atom ai */
1616 t[ai][0] = t[ai][0] + fix1;
1617 t[ai][1] = t[ai][1] + fiy1;
1618 t[ai][2] = t[ai][2] + fiz1;
1619
1620 fshift[shift][0] = fshift[shift][0] + fix1;
1621 fshift[shift][1] = fshift[shift][1] + fiy1;
1622 fshift[shift][2] = fshift[shift][2] + fiz1;
1623
1624 }
1625
1626 return 0;
1627 }
1628
1629
1630 real calc_gb_forces(t_commrec *cr, t_mdatoms *md, gmx_genborn_t *born, gmx_localtop_t *top, const t_atomtypes *atype,
1631 rvec x[], rvec f[], t_forcerec *fr, t_idef *idef, int gb_algorithm, t_nrnb *nrnb, bool bRad,
1632 const t_pbc *pbc, const t_graph *graph)
1633 {
1634 real v=0;
1635 int cnt;
1636
1637 /* PBC or not? */
1638 const t_pbc *pbc_null;
1639
1640 if (fr->bMolPBC)
1641 pbc_null = pbc;
1642 else
1643 pbc_null = NULL;
1644
1645
1646
1647 /* Do a simple ACE type approximation for the non-polar solvation */
1648 v += calc_gb_nonpolar(cr, fr,born->nr, born, top, atype, fr->dvda, gb_algorithm,md);
1649
1650 /* Calculate the bonded GB-interactions using either table or analytical formula */
1651 #ifdef GMX_DOUBLE
1652 v += gb_bonds_tab(x,f,fr->fshift, md->chargeA,&(fr->gbtabscale),
1653 fr->invsqrta,fr->dvda,fr->gbtab.tab,idef,born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph);
1654 #else
1655 #if ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || defined(GMX_SSE2) )
1656 v += gb_bonds_analytic(x[0],f[0],md->chargeA,born->bRad,fr->dvda,idef,born->gb_epsilon_solvent,fr->epsfac);
1657 #else
1658 v += gb_bonds_tab(x,f,fr->fshift, md->chargeA,&(fr->gbtabscale),
1659 fr->invsqrta,fr->dvda,fr->gbtab.tab,idef,born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph);
1660 #endif
1661 #endif
1662
1663 /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1664 v += calc_gb_selfcorrections(cr,born->nr,md->chargeA, born, fr->dvda, md, fr->epsfac);
1665
1666 /* If parallel, sum the derivative of the potential w.r.t the born radii */
1667 if(PARTDECOMP(cr))
1668 {
1669 gmx_sum(md->nr,fr->dvda, cr);
1670 }
1671 else if(DOMAINDECOMP(cr))
1672 {
1673 dd_atom_sum_real(cr->dd,fr->dvda);
1674 dd_atom_spread_real(cr->dd,fr->dvda);
1675 }
1676
1677 #ifndef DOUBLE
1678 if(fr->bAllvsAll)
1679 {
1680 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_SSE2) )
1681 genborn_allvsall_calc_chainrule_sse2_single(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1682 #else
1683 genborn_allvsall_calc_chainrule(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1684 #endif
1685 cnt = md->homenr*(md->nr/2+1);
1686 inc_nrnb(nrnb,eNR_BORN_AVA_CHAINRULE,cnt);
1687 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,md->homenr);
1688 return v;
1689 }
1690 #endif
1691
1692 #ifdef GMX_DOUBLE
1693
1694 #if ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || defined(GMX_SSE2) )
1695 calc_gb_chainrule_sse2_double(born->nr, &(fr->gblist), fr->dadx, fr->dvda,
1696 x[0], f[0], fr->fshift[0], fr->shift_vec[0],
1697 gb_algorithm, born);
1698 #else
1699 calc_gb_chainrule(born->nr, &(fr->gblist), fr->dadx, fr->dvda,
1700 x, f, fr->fshift, fr->shift_vec,
1701 gb_algorithm, born, md);
1702 #endif
1703
1704 #else
1705
1706 #if (!defined DISABLE_SSE && ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_SSE2) ))
1707 /* x86 or x86-64 with GCC inline assembly and/or SSE intrinsics */
1708 calc_gb_chainrule_sse(born->nr, &(fr->gblist), fr->dadx, fr->dvda,
1709 x[0], f[0], fr->fshift[0], fr->shift_vec[0],
1710 gb_algorithm, born, md);
1711 #else
1712 /* Calculate the forces due to chain rule terms with non sse code */
1713 calc_gb_chainrule(born->nr, &(fr->gblist), fr->dadx, fr->dvda,
1714 x, f, fr->fshift, fr->shift_vec,
1715 gb_algorithm, born, md);
1716 #endif
1717 #endif
1718
1719 if(!fr->bAllvsAll)
1720 {
1721 inc_nrnb(nrnb,eNR_BORN_CHAINRULE,fr->gblist.nrj);
1722 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,fr->gblist.nri);
1723
1724 }
1725
1726 return v;
1727
1728 }
1729
1730 static void add_j_to_gblist(gbtmpnbl_t *list,int aj)
1731 {
1732 if (list->naj >= list->aj_nalloc)
1733 {
1734 list->aj_nalloc = over_alloc_large(list->naj+1);
1735 srenew(list->aj,list->aj_nalloc);
1736 }
1737
1738 list->aj[list->naj++] = aj;
1739 }
1740
1741 static gbtmpnbl_t *find_gbtmplist(struct gbtmpnbls *lists,int shift)
1742 {
1743 int ind,i;
1744
1745 /* Search the list with the same shift, if there is one */
1746 ind = 0;
1747 while (ind < lists->nlist && shift != lists->list[ind].shift)
1748 {
1749 ind++;
1750 }
1751 if (ind == lists->nlist)
1752 {
1753 if (lists->nlist == lists->list_nalloc)
1754 {
1755 lists->list_nalloc++;
1756 srenew(lists->list,lists->list_nalloc);
1757 for(i=lists->nlist; i<lists->list_nalloc; i++)
1758 {
1759 lists->list[i].aj = NULL;
1760 lists->list[i].aj_nalloc = 0;
1761 }
1762
1763 }
1764
1765 lists->list[lists->nlist].shift = shift;
1766 lists->list[lists->nlist].naj = 0;
1767 lists->nlist++;
1768 }
1769
1770 return &lists->list[ind];
1771 }
1772
1773 static void add_bondeds_to_gblist(t_ilist *il,
1774 bool bMolPBC,t_pbc *pbc,t_graph *g,rvec *x,
1775 struct gbtmpnbls *nls)
1776 {
1777 int ind,j,ai,aj,shift,found;
1778 rvec dx;
1779 ivec dt;
1780 gbtmpnbl_t *list;
1781
1782 shift = CENTRAL;
1783 for(ind=0; ind<il->nr; ind+=3)
1784 {
1785 ai = il->iatoms[ind+1];
1786 aj = il->iatoms[ind+2];
1787
1788 shift = CENTRAL;
1789 if (g != NULL)
1790 {
1791 rvec_sub(x[ai],x[aj],dx);
1792 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1793 shift = IVEC2IS(dt);
1794 }
1795 else if (bMolPBC)
1796 {
1797 shift = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
1798 }
1799
1800 /* Find the list for this shift or create one */
1801 list = find_gbtmplist(&nls[ai],shift);
1802
1803 found=0;
1804
1805 /* So that we do not add the same bond twice.
1806 * This happens with some constraints between 1-3 atoms
1807 * that are in the bond-list but should not be in the GB nb-list */
1808 for(j=0;j<list->naj;j++)
1809 {
1810 if (list->aj[j] == aj)
1811 {
1812 found = 1;
1813 }
1814 }
1815
1816 if (found == 0)
1817 {
1818 if(ai == aj)
1819 {
1820 gmx_incons("ai == aj");
1821 }
1822
1823 add_j_to_gblist(list,aj);
1824 }
1825 }
1826 }
1827
1828 static int
1829 compare_int (const void * a, const void * b)
1830 {
1831 return ( *(int*)a - *(int*)b );
1832 }
1833
1834
1835
1836 int make_gb_nblist(t_commrec *cr, int gb_algorithm, real gbcut,
1837 rvec x[], matrix box,
1838 t_forcerec *fr, t_idef *idef, t_graph *graph, gmx_genborn_t *born)
1839 {
1840 int i,l,ii,j,k,n,nj0,nj1,ai,aj,at0,at1,found,shift,s;
1841 int apa;
1842 t_nblist *nblist;
1843 t_pbc pbc;
1844
1845 struct gbtmpnbls *nls;
1846 gbtmpnbl_t *list =NULL;
1847
1848 nls = born->nblist_work;
1849
1850 for(i=0;i<born->nr;i++)
1851 {
1852 nls[i].nlist = 0;
1853 }
1854
1855 if (fr->bMolPBC)
1856 {
1857 set_pbc_dd(&pbc,fr->ePBC,cr->dd,TRUE,box);
1858 }
1859
1860 switch (gb_algorithm)
1861 {
1862 case egbHCT:
1863 case egbOBC:
1864 /* Loop over 1-2, 1-3 and 1-4 interactions */
1865 for(j=F_GB12;j<=F_GB14;j++)
1866 {
1867 add_bondeds_to_gblist(&idef->il[j],fr->bMolPBC,&pbc,graph,x,nls);
1868 }
1869 break;
1870 case egbSTILL:
1871 /* Loop over 1-4 interactions */
1872 add_bondeds_to_gblist(&idef->il[F_GB14],fr->bMolPBC,&pbc,graph,x,nls);
1873 break;
1874 default:
1875 gmx_incons("Unknown GB algorithm");
1876 }
1877
1878 /* Loop over the VDWQQ and VDW nblists to set up the nonbonded part of the GB list */
1879 for(n=0; (n<fr->nnblists); n++)
1880 {
1881 for(i=0; (i<eNL_NR); i++)
1882 {
1883 nblist=&(fr->nblists[n].nlist_sr[i]);
1884
1885 if (nblist->nri > 0 && (i==eNL_VDWQQ || i==eNL_QQ))
1886 {
1887 for(j=0;j<nblist->nri;j++)
1888 {
1889 ai = nblist->iinr[j];
1890 shift = nblist->shift[j];
1891
1892 /* Find the list for this shift or create one */
1893 list = find_gbtmplist(&nls[ai],shift);
1894
1895 nj0 = nblist->jindex[j];
1896 nj1 = nblist->jindex[j+1];
1897
1898 /* Add all the j-atoms in the non-bonded list to the GB list */
1899 for(k=nj0;k<nj1;k++)
1900 {
1901 add_j_to_gblist(list,nblist->jjnr[k]);
1902 }
1903 }
1904 }
1905 }
1906 }
1907
1908 /* Zero out some counters */
1909 fr->gblist.nri=0;
1910 fr->gblist.nrj=0;
1911
1912 fr->gblist.jindex[0] = fr->gblist.nri;
1913
1914 for(i=0;i<fr->natoms_force;i++)
1915 {
1916 for(s=0; s<nls[i].nlist; s++)
1917 {
1918 list = &nls[i].list[s];
1919
1920 /* Only add those atoms that actually have neighbours */
1921 if (born->use[i] != 0)
1922 {
1923 fr->gblist.iinr[fr->gblist.nri] = i;
1924 fr->gblist.shift[fr->gblist.nri] = list->shift;
1925 fr->gblist.nri++;
1926
1927 for(k=0; k<list->naj; k++)
1928 {
1929 /* Memory allocation for jjnr */
1930 if(fr->gblist.nrj >= fr->gblist.maxnrj)
1931 {
1932 fr->gblist.maxnrj += over_alloc_large(fr->gblist.maxnrj);
1933
1934 if (debug)
1935 {
1936 fprintf(debug,"Increasing GB neighbourlist j size to %d\n",fr->gblist.maxnrj);
1937 }
1938
1939 srenew(fr->gblist.jjnr,fr->gblist.maxnrj);
1940 }
1941
1942 /* Put in list */
1943 if(i == list->aj[k])
1944 {
1945 gmx_incons("i == list->aj[k]");
1946 }
1947 fr->gblist.jjnr[fr->gblist.nrj++] = list->aj[k];
1948 }
1949
1950 fr->gblist.jindex[fr->gblist.nri] = fr->gblist.nrj;
1951 }
1952 }
1953 }
1954
1955
1956 #ifdef SORT_GB_LIST
1957 for(i=0;i<fr->gblist.nri;i++)
1958 {
1959 nj0 = fr->gblist.jindex[i];
1960 nj1 = fr->gblist.jindex[i+1];
1961 ai = fr->gblist.iinr[i];
1962
1963 /* Temporary fix */
1964 for(j=nj0;j<nj1;j++)
1965 {
1966 if(fr->gblist.jjnr[j]<ai)
1967 fr->gblist.jjnr[j]+=fr->natoms_force;
1968 }
1969 qsort(fr->gblist.jjnr+nj0,nj1-nj0,sizeof(int),compare_int);
1970 /* Fix back */
1971 for(j=nj0;j<nj1;j++)
1972 {
1973 if(fr->gblist.jjnr[j]>=fr->natoms_force)
1974 fr->gblist.jjnr[j]-=fr->natoms_force;
1975 }
1976
1977 }
1978 #endif
1979
1980 return 0;
1981 }
1982
1983 void make_local_gb(const t_commrec *cr, gmx_genborn_t *born, int gb_algorithm)
1984 {
1985 int i,at0,at1;
1986 gmx_domdec_t *dd=NULL;
1987
1988 if(DOMAINDECOMP(cr))
1989 {
1990 dd = cr->dd;
1991 at0 = 0;
1992 at1 = dd->nat_tot;
1993 }
1994 else
1995 {
1996 /* Single node or particle decomp (global==local), just copy pointers and return */
1997 if(gb_algorithm==egbSTILL)
1998 {
1999 born->gpol = born->gpol_globalindex;
2000 born->vsolv = born->vsolv_globalindex;
2001 born->gb_radius = born->gb_radius_globalindex;
2002 }
2003 else
2004 {
2005 born->param = born->param_globalindex;
2006 born->gb_radius = born->gb_radius_globalindex;
2007 }
2008
2009 born->use = born->use_globalindex;
2010
2011 return;
2012 }
2013
2014 /* Reallocation of local arrays if necessary */
2015 if(born->nlocal < dd->nat_tot)
2016 {
2017 born->nlocal = dd->nat_tot;
2018
2019 /* Arrays specific to different gb algorithms */
2020 if(gb_algorithm==egbSTILL)
2021 {
2022 srenew(born->gpol, born->nlocal+3);
2023 srenew(born->vsolv, born->nlocal+3);
2024 srenew(born->gb_radius, born->nlocal+3);
2025 }
2026 else
2027 {
2028 srenew(born->param, born->nlocal+3);
2029 srenew(born->gb_radius, born->nlocal+3);
2030 }
2031
2032 /* All gb-algorithms use the array for vsites exclusions */
2033 srenew(born->use, born->nlocal+3);
2034 }
2035
2036 /* With dd, copy algorithm specific arrays */
2037 if(gb_algorithm==egbSTILL)
2038 {
2039 for(i=at0;i<at1;i++)
2040 {
2041 born->gpol[i] = born->gpol_globalindex[dd->gatindex[i]];
2042 born->vsolv[i] = born->vsolv_globalindex[dd->gatindex[i]];
2043 born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
2044 born->use[i] = born->use_globalindex[dd->gatindex[i]];
2045 }
2046 }
2047 else
2048 {
2049 for(i=at0;i<at1;i++)
2050 {
2051 born->param[i] = born->param_globalindex[dd->gatindex[i]];
2052 born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
2053 born->use[i] = born->use_globalindex[dd->gatindex[i]];
2054 }
2055 }
2056 }
2057
Atempt to write QMMM interface for Gamess-US