3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2008, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
32 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
41 #include "gmx_fatal.h"
56 #include "mtop_util.h"
57 #include "chargegroup.h"
62 atom_id shell
; /* The shell id */
63 atom_id nucl1
,nucl2
,nucl3
; /* The nuclei connected to the shell */
64 /* gmx_bool bInterCG; */ /* Coupled to nuclei outside cg? */
65 real k
; /* force constant */
66 real k_1
; /* 1 over force constant */
72 typedef struct gmx_shellfc
{
73 int nshell_gl
; /* The number of shells in the system */
74 t_shell
*shell_gl
; /* All the shells (for DD only) */
75 int *shell_index_gl
; /* Global shell index (for DD only) */
76 gmx_bool bInterCG
; /* Are there inter charge-group shells? */
77 int nshell
; /* The number of local shells */
78 t_shell
*shell
; /* The local shells */
79 int shell_nalloc
; /* The allocation size of shell */
80 gmx_bool bPredict
; /* Predict shell positions */
81 gmx_bool bRequireInit
; /* Require initialization of shell positions */
82 int nflexcon
; /* The number of flexible constraints */
83 rvec
*x
[2]; /* Array for iterative minimization */
84 rvec
*f
[2]; /* Array for iterative minimization */
85 int x_nalloc
; /* The allocation size of x and f */
86 rvec
*acc_dir
; /* Acceleration direction for flexcon */
87 rvec
*x_old
; /* Old coordinates for flexcon */
88 int flex_nalloc
; /* The allocation size of acc_dir and x_old */
89 rvec
*adir_xnold
; /* Work space for init_adir */
90 rvec
*adir_xnew
; /* Work space for init_adir */
91 int adir_nalloc
; /* Work space for init_adir */
95 static void pr_shell(FILE *fplog
,int ns
,t_shell s
[])
99 fprintf(fplog
,"SHELL DATA\n");
100 fprintf(fplog
,"%5s %8s %5s %5s %5s\n",
101 "Shell","Force k","Nucl1","Nucl2","Nucl3");
102 for(i
=0; (i
<ns
); i
++) {
103 fprintf(fplog
,"%5d %8.3f %5d",s
[i
].shell
,1.0/s
[i
].k_1
,s
[i
].nucl1
);
105 fprintf(fplog
," %5d\n",s
[i
].nucl2
);
106 else if (s
[i
].nnucl
== 3)
107 fprintf(fplog
," %5d %5d\n",s
[i
].nucl2
,s
[i
].nucl3
);
113 static void predict_shells(FILE *fplog
,rvec x
[],rvec v
[],real dt
,
115 real mass
[],gmx_mtop_t
*mtop
,gmx_bool bInit
)
118 real dt_1
,dt_2
,dt_3
,fudge
,tm
,m1
,m2
,m3
;
120 gmx_mtop_atomlookup_t alook
=NULL
;
124 alook
= gmx_mtop_atomlookup_init(mtop
);
127 /* We introduce a fudge factor for performance reasons: with this choice
128 * the initial force on the shells is about a factor of two lower than
135 fprintf(fplog
,"RELAX: Using prediction for initial shell placement\n");
144 for(i
=0; (i
<ns
); i
++) {
148 switch (s
[i
].nnucl
) {
151 for(m
=0; (m
<DIM
); m
++)
152 x
[s1
][m
]+=ptr
[n1
][m
]*dt_1
;
161 /* Not the correct masses with FE, but it is just a prediction... */
166 for(m
=0; (m
<DIM
); m
++)
167 x
[s1
][m
]+=(m1
*ptr
[n1
][m
]+m2
*ptr
[n2
][m
])*tm
;
178 /* Not the correct masses with FE, but it is just a prediction... */
179 gmx_mtop_atomnr_to_atom(alook
,n1
,&atom
);
181 gmx_mtop_atomnr_to_atom(alook
,n2
,&atom
);
183 gmx_mtop_atomnr_to_atom(alook
,n3
,&atom
);
186 tm
= dt_1
/(m1
+m2
+m3
);
187 for(m
=0; (m
<DIM
); m
++)
188 x
[s1
][m
]+=(m1
*ptr
[n1
][m
]+m2
*ptr
[n2
][m
]+m3
*ptr
[n3
][m
])*tm
;
191 gmx_fatal(FARGS
,"Shell %d has %d nuclei!",i
,s
[i
].nnucl
);
196 gmx_mtop_atomlookup_destroy(alook
);
200 gmx_shellfc_t
init_shell_flexcon(FILE *fplog
,
201 gmx_mtop_t
*mtop
,int nflexcon
,
204 struct gmx_shellfc
*shfc
;
206 int *shell_index
=NULL
,*at2cg
;
208 int n
[eptNR
],ns
,nshell
,nsi
;
209 int i
,j
,nmol
,type
,mb
,mt
,a_offset
,cg
,mol
,ftype
,nra
;
211 int aS
,aN
=0; /* Shell and nucleus */
212 int bondtypes
[] = { F_BONDS
, F_HARMONIC
, F_CUBICBONDS
, F_POLARIZATION
, F_ANHARM_POL
, F_WATER_POL
};
213 #define NBT asize(bondtypes)
215 gmx_mtop_atomloop_block_t aloopb
;
216 gmx_mtop_atomloop_all_t aloop
;
217 gmx_ffparams_t
*ffparams
;
218 gmx_molblock_t
*molb
;
222 /* Count number of shells, and find their indices */
223 for(i
=0; (i
<eptNR
); i
++) {
227 aloopb
= gmx_mtop_atomloop_block_init(mtop
);
228 while (gmx_mtop_atomloop_block_next(aloopb
,&atom
,&nmol
)) {
229 n
[atom
->ptype
] += nmol
;
233 /* Print the number of each particle type */
234 for(i
=0; (i
<eptNR
); i
++) {
236 fprintf(fplog
,"There are: %d %ss\n",n
[i
],ptype_str
[i
]);
241 nshell
= n
[eptShell
];
243 if (nshell
== 0 && nflexcon
== 0) {
248 shfc
->nflexcon
= nflexcon
;
254 /* We have shells: fill the shell data structure */
256 /* Global system sized array, this should be avoided */
257 snew(shell_index
,mtop
->natoms
);
259 aloop
= gmx_mtop_atomloop_all_init(mtop
);
261 while (gmx_mtop_atomloop_all_next(aloop
,&i
,&atom
)) {
262 if (atom
->ptype
== eptShell
) {
263 shell_index
[i
] = nshell
++;
269 /* Initiate the shell structures */
270 for(i
=0; (i
<nshell
); i
++) {
271 shell
[i
].shell
= NO_ATID
;
273 shell
[i
].nucl1
= NO_ATID
;
274 shell
[i
].nucl2
= NO_ATID
;
275 shell
[i
].nucl3
= NO_ATID
;
276 /* shell[i].bInterCG=FALSE; */
281 ffparams
= &mtop
->ffparams
;
283 /* Now fill the structures */
284 shfc
->bInterCG
= FALSE
;
287 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
288 molb
= &mtop
->molblock
[mb
];
289 molt
= &mtop
->moltype
[molb
->type
];
292 snew(at2cg
,molt
->atoms
.nr
);
293 for(cg
=0; cg
<cgs
->nr
; cg
++) {
294 for(i
=cgs
->index
[cg
]; i
<cgs
->index
[cg
+1]; i
++) {
299 atom
= molt
->atoms
.atom
;
300 for(mol
=0; mol
<molb
->nmol
; mol
++) {
301 for(j
=0; (j
<NBT
); j
++) {
302 ia
= molt
->ilist
[bondtypes
[j
]].iatoms
;
303 for(i
=0; (i
<molt
->ilist
[bondtypes
[j
]].nr
); ) {
305 ftype
= ffparams
->functype
[type
];
306 nra
= interaction_function
[ftype
].nratoms
;
308 /* Check whether we have a bond with a shell */
311 switch (bondtypes
[j
]) {
317 if (atom
[ia
[1]].ptype
== eptShell
) {
321 else if (atom
[ia
[2]].ptype
== eptShell
) {
327 aN
= ia
[4]; /* Dummy */
328 aS
= ia
[5]; /* Shell */
331 gmx_fatal(FARGS
,"Death Horror: %s, %d",__FILE__
,__LINE__
);
337 /* Check whether one of the particles is a shell... */
338 nsi
= shell_index
[a_offset
+aS
];
339 if ((nsi
< 0) || (nsi
>= nshell
))
340 gmx_fatal(FARGS
,"nsi is %d should be within 0 - %d. aS = %d",
342 if (shell
[nsi
].shell
== NO_ATID
) {
343 shell
[nsi
].shell
= a_offset
+ aS
;
346 else if (shell
[nsi
].shell
!= a_offset
+aS
)
347 gmx_fatal(FARGS
,"Weird stuff in %s, %d",__FILE__
,__LINE__
);
349 if (shell
[nsi
].nucl1
== NO_ATID
) {
350 shell
[nsi
].nucl1
= a_offset
+ aN
;
351 } else if (shell
[nsi
].nucl2
== NO_ATID
) {
352 shell
[nsi
].nucl2
= a_offset
+ aN
;
353 } else if (shell
[nsi
].nucl3
== NO_ATID
) {
354 shell
[nsi
].nucl3
= a_offset
+ aN
;
357 pr_shell(fplog
,ns
,shell
);
358 gmx_fatal(FARGS
,"Can not handle more than three bonds per shell\n");
360 if (at2cg
[aS
] != at2cg
[aN
]) {
361 /* shell[nsi].bInterCG = TRUE; */
362 shfc
->bInterCG
= TRUE
;
365 switch (bondtypes
[j
]) {
368 shell
[nsi
].k
+= ffparams
->iparams
[type
].harmonic
.krA
;
371 shell
[nsi
].k
+= ffparams
->iparams
[type
].cubic
.kb
;
375 if (!gmx_within_tol(qS
, atom
[aS
].qB
, GMX_REAL_EPS
*10))
376 gmx_fatal(FARGS
,"polarize can not be used with qA(%e) != qB(%e) for atom %d of molecule block %d", qS
, atom
[aS
].qB
, aS
+1, mb
+1);
377 shell
[nsi
].k
+= sqr(qS
)*ONE_4PI_EPS0
/
378 ffparams
->iparams
[type
].polarize
.alpha
;
381 if (!gmx_within_tol(qS
, atom
[aS
].qB
, GMX_REAL_EPS
*10))
382 gmx_fatal(FARGS
,"water_pol can not be used with qA(%e) != qB(%e) for atom %d of molecule block %d", qS
, atom
[aS
].qB
, aS
+1, mb
+1);
383 alpha
= (ffparams
->iparams
[type
].wpol
.al_x
+
384 ffparams
->iparams
[type
].wpol
.al_y
+
385 ffparams
->iparams
[type
].wpol
.al_z
)/3.0;
386 shell
[nsi
].k
+= sqr(qS
)*ONE_4PI_EPS0
/alpha
;
389 gmx_fatal(FARGS
,"Death Horror: %s, %d",__FILE__
,__LINE__
);
397 a_offset
+= molt
->atoms
.nr
;
399 /* Done with this molecule type */
403 /* Verify whether it's all correct */
405 gmx_fatal(FARGS
,"Something weird with shells. They may not be bonded to something");
407 for(i
=0; (i
<ns
); i
++)
408 shell
[i
].k_1
= 1.0/shell
[i
].k
;
411 pr_shell(debug
,ns
,shell
);
414 shfc
->nshell_gl
= ns
;
415 shfc
->shell_gl
= shell
;
416 shfc
->shell_index_gl
= shell_index
;
418 shfc
->bPredict
= (getenv("GMX_NOPREDICT") == NULL
);
419 shfc
->bRequireInit
= FALSE
;
420 if (!shfc
->bPredict
) {
422 fprintf(fplog
,"\nWill never predict shell positions\n");
424 shfc
->bRequireInit
= (getenv("GMX_REQUIRE_SHELL_INIT") != NULL
);
425 if (shfc
->bRequireInit
&& fplog
)
426 fprintf(fplog
,"\nWill always initiate shell positions\n");
429 if (shfc
->bPredict
) {
431 predict_shells(fplog
,x
,NULL
,0,shfc
->nshell_gl
,shfc
->shell_gl
,
435 if (shfc
->bInterCG
) {
437 fprintf(fplog
,"\nNOTE: there all shells that are connected to particles outside thier own charge group, will not predict shells positions during the run\n\n");
438 shfc
->bPredict
= FALSE
;
445 void make_local_shells(t_commrec
*cr
,t_mdatoms
*md
,
446 struct gmx_shellfc
*shfc
)
449 int a0
,a1
,*ind
,nshell
,i
;
450 gmx_domdec_t
*dd
=NULL
;
453 if (DOMAINDECOMP(cr
)) {
458 pd_at_range(cr
,&a0
,&a1
);
461 /* Single node: we need all shells, just copy the pointer */
462 shfc
->nshell
= shfc
->nshell_gl
;
463 shfc
->shell
= shfc
->shell_gl
;
468 ind
= shfc
->shell_index_gl
;
472 for(i
=a0
; i
<a1
; i
++) {
473 if (md
->ptype
[i
] == eptShell
) {
474 if (nshell
+1 > shfc
->shell_nalloc
) {
475 shfc
->shell_nalloc
= over_alloc_dd(nshell
+1);
476 srenew(shell
,shfc
->shell_nalloc
);
479 shell
[nshell
] = shfc
->shell_gl
[ind
[dd
->gatindex
[i
]]];
481 shell
[nshell
] = shfc
->shell_gl
[ind
[i
]];
483 /* With inter-cg shells we can no do shell prediction,
484 * so we do not need the nuclei numbers.
486 if (!shfc
->bInterCG
) {
487 shell
[nshell
].nucl1
= i
+ shell
[nshell
].nucl1
- shell
[nshell
].shell
;
488 if (shell
[nshell
].nnucl
> 1)
489 shell
[nshell
].nucl2
= i
+ shell
[nshell
].nucl2
- shell
[nshell
].shell
;
490 if (shell
[nshell
].nnucl
> 2)
491 shell
[nshell
].nucl3
= i
+ shell
[nshell
].nucl3
- shell
[nshell
].shell
;
493 shell
[nshell
].shell
= i
;
498 shfc
->nshell
= nshell
;
502 static void do_1pos(rvec xnew
,rvec xold
,rvec f
,real step
)
520 static void do_1pos3(rvec xnew
,rvec xold
,rvec f
,rvec step
)
538 static void directional_sd(FILE *log
,rvec xold
[],rvec xnew
[],rvec acc_dir
[],
539 int start
,int homenr
,real step
)
543 for(i
=start
; i
<homenr
; i
++)
544 do_1pos(xnew
[i
],xold
[i
],acc_dir
[i
],step
);
547 static void shell_pos_sd(FILE *log
,rvec xcur
[],rvec xnew
[],rvec f
[],
548 int ns
,t_shell s
[],int count
)
550 const real step_scale_min
= 0.8,
551 step_scale_increment
= 0.2,
552 step_scale_max
= 1.2,
553 step_scale_multiple
= (step_scale_max
- step_scale_min
) / step_scale_increment
;
557 real step_min
,step_max
;
562 for(i
=0; (i
<ns
); i
++) {
565 for(d
=0; d
<DIM
; d
++) {
566 s
[i
].step
[d
] = s
[i
].k_1
;
568 step_min
= min(step_min
,s
[i
].step
[d
]);
569 step_max
= max(step_max
,s
[i
].step
[d
]);
573 for(d
=0; d
<DIM
; d
++) {
574 dx
= xcur
[shell
][d
] - s
[i
].xold
[d
];
575 df
= f
[shell
][d
] - s
[i
].fold
[d
];
576 /* -dx/df gets used to generate an interpolated value, but would
577 * cause a NaN if df were binary-equal to zero. Values close to
578 * zero won't cause problems (because of the min() and max()), so
579 * just testing for binary inequality is OK. */
583 /* Scale the step size by a factor interpolated from
584 * step_scale_min to step_scale_max, as k_est goes from 0 to
585 * step_scale_multiple * s[i].step[d] */
587 step_scale_min
* s
[i
].step
[d
] +
588 step_scale_increment
* min(step_scale_multiple
* s
[i
].step
[d
], max(k_est
, 0));
593 if (gmx_numzero(dx
)) /* 0 == dx */
595 /* Likely this will never happen, but if it does just
596 * don't scale the step. */
600 s
[i
].step
[d
] *= step_scale_max
;
604 step_min
= min(step_min
,s
[i
].step
[d
]);
605 step_max
= max(step_max
,s
[i
].step
[d
]);
609 copy_rvec(xcur
[shell
],s
[i
].xold
);
610 copy_rvec(f
[shell
], s
[i
].fold
);
612 do_1pos3(xnew
[shell
],xcur
[shell
],f
[shell
],s
[i
].step
);
615 fprintf(debug
,"shell[%d] = %d\n",i
,shell
);
616 pr_rvec(debug
,0,"fshell",f
[shell
],DIM
,TRUE
);
617 pr_rvec(debug
,0,"xold",xcur
[shell
],DIM
,TRUE
);
618 pr_rvec(debug
,0,"step",s
[i
].step
,DIM
,TRUE
);
619 pr_rvec(debug
,0,"xnew",xnew
[shell
],DIM
,TRUE
);
623 printf("step %.3e %.3e\n",step_min
,step_max
);
627 static void decrease_step_size(int nshell
,t_shell s
[])
631 for(i
=0; i
<nshell
; i
++)
632 svmul(0.8,s
[i
].step
,s
[i
].step
);
635 static void print_epot(FILE *fp
,gmx_large_int_t mdstep
,int count
,real epot
,real df
,
636 int ndir
,real sf_dir
)
640 fprintf(fp
,"MDStep=%5s/%2d EPot: %12.8e, rmsF: %6.2e",
641 gmx_step_str(mdstep
,buf
),count
,epot
,df
);
643 fprintf(fp
,", dir. rmsF: %6.2e\n",sqrt(sf_dir
/ndir
));
649 static real
rms_force(t_commrec
*cr
,rvec f
[],int ns
,t_shell s
[],
650 int ndir
,real
*sf_dir
,real
*Epot
)
656 for(i
=0; i
<ns
; i
++) {
658 buf
[0] += norm2(f
[shell
]);
667 ntot
= (int)(buf
[1] + 0.5);
673 return (ntot
? sqrt(buf
[0]/ntot
) : 0);
676 static void check_pbc(FILE *fp
,rvec x
[],int shell
)
681 for(m
=0; (m
<DIM
); m
++)
682 if (fabs(x
[shell
][m
]-x
[now
][m
]) > 0.3) {
683 pr_rvecs(fp
,0,"SHELL-X",x
+now
,5);
688 static void dump_shells(FILE *fp
,rvec x
[],rvec f
[],real ftol
,int ns
,t_shell s
[])
695 for(i
=0; (i
<ns
); i
++) {
697 ff2
= iprod(f
[shell
],f
[shell
]);
699 fprintf(fp
,"SHELL %5d, force %10.5f %10.5f %10.5f, |f| %10.5f\n",
700 shell
,f
[shell
][XX
],f
[shell
][YY
],f
[shell
][ZZ
],sqrt(ff2
));
701 check_pbc(fp
,x
,shell
);
705 static void init_adir(FILE *log
,gmx_shellfc_t shfc
,
706 gmx_constr_t constr
,t_idef
*idef
,t_inputrec
*ir
,
707 t_commrec
*cr
,int dd_ac1
,
708 gmx_large_int_t step
,t_mdatoms
*md
,int start
,int end
,
709 rvec
*x_old
,rvec
*x_init
,rvec
*x
,
710 rvec
*f
,rvec
*acc_dir
,
711 gmx_bool bMolPBC
,matrix box
,
712 real
*lambda
,real
*dvdlambda
,t_nrnb
*nrnb
)
719 unsigned short *ptype
;
722 if (DOMAINDECOMP(cr
))
726 if (n
> shfc
->adir_nalloc
) {
727 shfc
->adir_nalloc
= over_alloc_dd(n
);
728 srenew(shfc
->adir_xnold
,shfc
->adir_nalloc
);
729 srenew(shfc
->adir_xnew
,shfc
->adir_nalloc
);
731 xnold
= shfc
->adir_xnold
;
732 xnew
= shfc
->adir_xnew
;
738 /* Does NOT work with freeze or acceleration groups (yet) */
739 for (n
=start
; n
<end
; n
++) {
740 w_dt
= md
->invmass
[n
]*dt
;
742 for (d
=0; d
<DIM
; d
++) {
743 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
)) {
744 xnold
[n
-start
][d
] = x
[n
][d
] - (x_init
[n
][d
] - x_old
[n
][d
]);
745 xnew
[n
-start
][d
] = 2*x
[n
][d
] - x_old
[n
][d
] + f
[n
][d
]*w_dt
*dt
;
747 xnold
[n
-start
][d
] = x
[n
][d
];
748 xnew
[n
-start
][d
] = x
[n
][d
];
752 constrain(log
,FALSE
,FALSE
,constr
,idef
,ir
,NULL
,cr
,step
,0,md
,
753 x
,xnold
-start
,NULL
,bMolPBC
,box
,
754 lambda
[efptBONDED
],&(dvdlambda
[efptBONDED
]),
755 NULL
,NULL
,nrnb
,econqCoord
,FALSE
,0,0);
756 constrain(log
,FALSE
,FALSE
,constr
,idef
,ir
,NULL
,cr
,step
,0,md
,
757 x
,xnew
-start
,NULL
,bMolPBC
,box
,
758 lambda
[efptBONDED
],&(dvdlambda
[efptBONDED
]),
759 NULL
,NULL
,nrnb
,econqCoord
,FALSE
,0,0);
761 for (n
=start
; n
<end
; n
++) {
764 -(2*x
[n
][d
]-xnold
[n
-start
][d
]-xnew
[n
-start
][d
])/sqr(dt
)
765 - f
[n
][d
]*md
->invmass
[n
];
766 clear_rvec(acc_dir
[n
]);
769 /* Project the acceleration on the old bond directions */
770 constrain(log
,FALSE
,FALSE
,constr
,idef
,ir
,NULL
,cr
,step
,0,md
,
771 x_old
,xnew
-start
,acc_dir
,bMolPBC
,box
,
772 lambda
[efptBONDED
],&(dvdlambda
[efptBONDED
]),
773 NULL
,NULL
,nrnb
,econqDeriv_FlexCon
,FALSE
,0,0);
776 int relax_shell_flexcon(FILE *fplog
,t_commrec
*cr
,gmx_bool bVerbose
,
777 gmx_large_int_t mdstep
,t_inputrec
*inputrec
,
778 gmx_bool bDoNS
,int force_flags
,
783 gmx_enerdata_t
*enerd
,t_fcdata
*fcd
,
784 t_state
*state
,rvec f
[],
787 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
789 gmx_groups_t
*groups
,
790 struct gmx_shellfc
*shfc
,
793 double t
,rvec mu_tot
,
794 int natoms
,gmx_bool
*bConverged
,
801 rvec
*pos
[2],*force
[2],*acc_dir
=NULL
,*x_old
=NULL
;
805 real ftol
,xiH
,xiS
,dum
=0;
807 gmx_bool bCont
,bInit
;
808 int nat
,dd_ac0
,dd_ac1
=0,i
;
809 int start
=md
->start
,homenr
=md
->homenr
,end
=start
+homenr
,cg0
,cg1
;
810 int nflexcon
,g
,number_steps
,d
,Min
=0,count
=0;
811 #define Try (1-Min) /* At start Try = 1 */
813 bCont
= (mdstep
== inputrec
->init_step
) && inputrec
->bContinuation
;
814 bInit
= (mdstep
== inputrec
->init_step
) || shfc
->bRequireInit
;
815 ftol
= inputrec
->em_tol
;
816 number_steps
= inputrec
->niter
;
817 nshell
= shfc
->nshell
;
819 nflexcon
= shfc
->nflexcon
;
823 if (DOMAINDECOMP(cr
)) {
824 nat
= dd_natoms_vsite(cr
->dd
);
826 dd_get_constraint_range(cr
->dd
,&dd_ac0
,&dd_ac1
);
827 nat
= max(nat
,dd_ac1
);
833 if (nat
> shfc
->x_nalloc
) {
834 /* Allocate local arrays */
835 shfc
->x_nalloc
= over_alloc_dd(nat
);
836 for(i
=0; (i
<2); i
++) {
837 srenew(shfc
->x
[i
],shfc
->x_nalloc
);
838 srenew(shfc
->f
[i
],shfc
->x_nalloc
);
841 for(i
=0; (i
<2); i
++) {
843 force
[i
] = shfc
->f
[i
];
846 /* With particle decomposition this code only works
847 * when all particles involved with each shell are in the same cg.
850 if (bDoNS
&& inputrec
->ePBC
!= epbcNONE
&& !DOMAINDECOMP(cr
)) {
851 /* This is the only time where the coordinates are used
852 * before do_force is called, which normally puts all
853 * charge groups in the box.
855 if (PARTDECOMP(cr
)) {
856 pd_cg_range(cr
,&cg0
,&cg1
);
861 put_charge_groups_in_box(fplog
,cg0
,cg1
,fr
->ePBC
,state
->box
,
862 &(top
->cgs
),state
->x
,fr
->cg_cm
);
864 mk_mshift(fplog
,graph
,fr
->ePBC
,state
->box
,state
->x
);
867 /* After this all coordinate arrays will contain whole molecules */
869 shift_self(graph
,state
->box
,state
->x
);
872 if (nat
> shfc
->flex_nalloc
) {
873 shfc
->flex_nalloc
= over_alloc_dd(nat
);
874 srenew(shfc
->acc_dir
,shfc
->flex_nalloc
);
875 srenew(shfc
->x_old
,shfc
->flex_nalloc
);
877 acc_dir
= shfc
->acc_dir
;
879 for(i
=0; i
<homenr
; i
++) {
882 state
->x
[start
+i
][d
] - state
->v
[start
+i
][d
]*inputrec
->delta_t
;
886 /* Do a prediction of the shell positions */
887 if (shfc
->bPredict
&& !bCont
) {
888 predict_shells(fplog
,state
->x
,state
->v
,inputrec
->delta_t
,nshell
,shell
,
889 md
->massT
,NULL
,bInit
);
892 /* do_force expected the charge groups to be in the box */
894 unshift_self(graph
,state
->box
,state
->x
);
896 /* Calculate the forces first time around */
898 pr_rvecs(debug
,0,"x b4 do_force",state
->x
+ start
,homenr
);
900 do_force(fplog
,cr
,inputrec
,mdstep
,nrnb
,wcycle
,top
,mtop
,groups
,
901 state
->box
,state
->x
,&state
->hist
,
902 force
[Min
],force_vir
,md
,enerd
,fcd
,
904 fr
,vsite
,mu_tot
,t
,fp_field
,NULL
,bBornRadii
,
905 (bDoNS
? GMX_FORCE_NS
: 0) | force_flags
);
909 init_adir(fplog
,shfc
,
910 constr
,idef
,inputrec
,cr
,dd_ac1
,mdstep
,md
,start
,end
,
911 shfc
->x_old
-start
,state
->x
,state
->x
,force
[Min
],
913 fr
->bMolPBC
,state
->box
,state
->lambda
,&dum
,nrnb
);
915 for(i
=start
; i
<end
; i
++)
916 sf_dir
+= md
->massT
[i
]*norm2(shfc
->acc_dir
[i
-start
]);
919 Epot
[Min
] = enerd
->term
[F_EPOT
];
921 df
[Min
]=rms_force(cr
,shfc
->f
[Min
],nshell
,shell
,nflexcon
,&sf_dir
,&Epot
[Min
]);
924 fprintf(debug
,"df = %g %g\n",df
[Min
],df
[Try
]);
928 pr_rvecs(debug
,0,"force0",force
[Min
],md
->nr
);
931 if (nshell
+nflexcon
> 0) {
932 /* Copy x to pos[Min] & pos[Try]: during minimization only the
933 * shell positions are updated, therefore the other particles must
936 memcpy(pos
[Min
],state
->x
,nat
*sizeof(state
->x
[0]));
937 memcpy(pos
[Try
],state
->x
,nat
*sizeof(state
->x
[0]));
940 if (bVerbose
&& MASTER(cr
))
941 print_epot(stdout
,mdstep
,0,Epot
[Min
],df
[Min
],nflexcon
,sf_dir
);
944 fprintf(debug
,"%17s: %14.10e\n",
945 interaction_function
[F_EKIN
].longname
,enerd
->term
[F_EKIN
]);
946 fprintf(debug
,"%17s: %14.10e\n",
947 interaction_function
[F_EPOT
].longname
,enerd
->term
[F_EPOT
]);
948 fprintf(debug
,"%17s: %14.10e\n",
949 interaction_function
[F_ETOT
].longname
,enerd
->term
[F_ETOT
]);
950 fprintf(debug
,"SHELLSTEP %s\n",gmx_step_str(mdstep
,sbuf
));
953 /* First check whether we should do shells, or whether the force is
954 * low enough even without minimization.
956 *bConverged
= (df
[Min
] < ftol
);
958 for(count
=1; (!(*bConverged
) && (count
< number_steps
)); count
++) {
960 construct_vsites(fplog
,vsite
,pos
[Min
],nrnb
,inputrec
->delta_t
,state
->v
,
961 idef
->iparams
,idef
->il
,
962 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
965 init_adir(fplog
,shfc
,
966 constr
,idef
,inputrec
,cr
,dd_ac1
,mdstep
,md
,start
,end
,
967 x_old
-start
,state
->x
,pos
[Min
],force
[Min
],acc_dir
-start
,
968 fr
->bMolPBC
,state
->box
,state
->lambda
,&dum
,nrnb
);
970 directional_sd(fplog
,pos
[Min
],pos
[Try
],acc_dir
-start
,start
,end
,
974 /* New positions, Steepest descent */
975 shell_pos_sd(fplog
,pos
[Min
],pos
[Try
],force
[Min
],nshell
,shell
,count
);
977 /* do_force expected the charge groups to be in the box */
979 unshift_self(graph
,state
->box
,pos
[Try
]);
982 pr_rvecs(debug
,0,"RELAX: pos[Min] ",pos
[Min
] + start
,homenr
);
983 pr_rvecs(debug
,0,"RELAX: pos[Try] ",pos
[Try
] + start
,homenr
);
985 /* Try the new positions */
986 do_force(fplog
,cr
,inputrec
,1,nrnb
,wcycle
,
987 top
,mtop
,groups
,state
->box
,pos
[Try
],&state
->hist
,
988 force
[Try
],force_vir
,
989 md
,enerd
,fcd
,state
->lambda
,graph
,
990 fr
,vsite
,mu_tot
,t
,fp_field
,NULL
,bBornRadii
,
994 pr_rvecs(debug
,0,"RELAX: force[Min]",force
[Min
] + start
,homenr
);
995 pr_rvecs(debug
,0,"RELAX: force[Try]",force
[Try
] + start
,homenr
);
999 init_adir(fplog
,shfc
,
1000 constr
,idef
,inputrec
,cr
,dd_ac1
,mdstep
,md
,start
,end
,
1001 x_old
-start
,state
->x
,pos
[Try
],force
[Try
],acc_dir
-start
,
1002 fr
->bMolPBC
,state
->box
,state
->lambda
,&dum
,nrnb
);
1004 for(i
=start
; i
<end
; i
++)
1005 sf_dir
+= md
->massT
[i
]*norm2(acc_dir
[i
-start
]);
1008 Epot
[Try
] = enerd
->term
[F_EPOT
];
1010 df
[Try
]=rms_force(cr
,force
[Try
],nshell
,shell
,nflexcon
,&sf_dir
,&Epot
[Try
]);
1013 fprintf(debug
,"df = %g %g\n",df
[Min
],df
[Try
]);
1017 pr_rvecs(debug
,0,"F na do_force",force
[Try
] + start
,homenr
);
1019 fprintf(debug
,"SHELL ITER %d\n",count
);
1020 dump_shells(debug
,pos
[Try
],force
[Try
],ftol
,nshell
,shell
);
1024 if (bVerbose
&& MASTER(cr
))
1025 print_epot(stdout
,mdstep
,count
,Epot
[Try
],df
[Try
],nflexcon
,sf_dir
);
1027 *bConverged
= (df
[Try
] < ftol
);
1029 if ((df
[Try
] < df
[Min
])) {
1031 fprintf(debug
,"Swapping Min and Try\n");
1033 /* Correct the velocities for the flexible constraints */
1034 invdt
= 1/inputrec
->delta_t
;
1035 for(i
=start
; i
<end
; i
++) {
1036 for(d
=0; d
<DIM
; d
++)
1037 state
->v
[i
][d
] += (pos
[Try
][i
][d
] - pos
[Min
][i
][d
])*invdt
;
1042 decrease_step_size(nshell
,shell
);
1045 if (MASTER(cr
) && !(*bConverged
)) {
1046 /* Note that the energies and virial are incorrect when not converged */
1049 "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1050 gmx_step_str(mdstep
,sbuf
),number_steps
,df
[Min
]);
1052 "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1053 gmx_step_str(mdstep
,sbuf
),number_steps
,df
[Min
]);
1056 /* Copy back the coordinates and the forces */
1057 memcpy(state
->x
,pos
[Min
],nat
*sizeof(state
->x
[0]));
1058 memcpy(f
,force
[Min
],nat
*sizeof(f
[0]));