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