641ed8e34206d3448beb84cbee3daa0c2d257968
[alexxy/gromacs.git] / src / gromacs / gmxlib / typedefs.c
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5  * Copyright (c) 2001-2004, The GROMACS development team.
6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8  * and including many others, as listed in the AUTHORS file in the
9  * top-level source directory and at http://www.gromacs.org.
10  *
11  * GROMACS is free software; you can redistribute it and/or
12  * modify it under the terms of the GNU Lesser General Public License
13  * as published by the Free Software Foundation; either version 2.1
14  * of the License, or (at your option) any later version.
15  *
16  * GROMACS is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19  * Lesser General Public License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with GROMACS; if not, see
23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
25  *
26  * If you want to redistribute modifications to GROMACS, please
27  * consider that scientific software is very special. Version
28  * control is crucial - bugs must be traceable. We will be happy to
29  * consider code for inclusion in the official distribution, but
30  * derived work must not be called official GROMACS. Details are found
31  * in the README & COPYING files - if they are missing, get the
32  * official version at http://www.gromacs.org.
33  *
34  * To help us fund GROMACS development, we humbly ask that you cite
35  * the research papers on the package. Check out http://www.gromacs.org.
36  */
37 /* This file is completely threadsafe - keep it that way! */
38 #include "typedefs.h"
39
40 #ifdef HAVE_CONFIG_H
41 #include <config.h>
42 #endif
43
44 #include <string.h>
45
46 #include "thread_mpi/threads.h"
47
48 #include "macros.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/pbcutil/pbc.h"
51 #include "gromacs/random/random.h"
52 #include "gromacs/topology/atoms.h"
53 #include "gromacs/topology/block.h"
54 #include "gromacs/topology/symtab.h"
55 #include "gromacs/utility/smalloc.h"
56
57 /* The source code in this file should be thread-safe.
58       Please keep it that way. */
59
60
61
62 static gmx_bool            bOverAllocDD     = FALSE;
63 static tMPI_Thread_mutex_t over_alloc_mutex = TMPI_THREAD_MUTEX_INITIALIZER;
64
65
66 void set_over_alloc_dd(gmx_bool set)
67 {
68     tMPI_Thread_mutex_lock(&over_alloc_mutex);
69     /* we just make sure that we don't set this at the same time.
70        We don't worry too much about reading this rarely-set variable */
71     bOverAllocDD = set;
72     tMPI_Thread_mutex_unlock(&over_alloc_mutex);
73 }
74
75 int over_alloc_dd(int n)
76 {
77     if (bOverAllocDD)
78     {
79         return OVER_ALLOC_FAC*n + 100;
80     }
81     else
82     {
83         return n;
84     }
85 }
86
87 int gmx_int64_to_int(gmx_int64_t step, const char *warn)
88 {
89     int i;
90
91     i = (int)step;
92
93     if (warn != NULL && (step < INT_MIN || step > INT_MAX))
94     {
95         fprintf(stderr, "\nWARNING during %s:\n", warn);
96         fprintf(stderr, "step value ");
97         fprintf(stderr, "%"GMX_PRId64, step);
98         fprintf(stderr, " does not fit in int, converted to %d\n\n", i);
99     }
100
101     return i;
102 }
103
104 char *gmx_step_str(gmx_int64_t i, char *buf)
105 {
106     sprintf(buf, "%"GMX_PRId64, i);
107
108     return buf;
109 }
110
111 void init_groups(gmx_groups_t *groups)
112 {
113     int g;
114
115     groups->ngrpname = 0;
116     groups->grpname  = NULL;
117     for (g = 0; (g < egcNR); g++)
118     {
119         groups->grps[g].nm_ind = NULL;
120         groups->ngrpnr[g]      = 0;
121         groups->grpnr[g]       = NULL;
122     }
123
124 }
125
126 void init_mtop(gmx_mtop_t *mtop)
127 {
128     mtop->name         = NULL;
129     mtop->nmoltype     = 0;
130     mtop->moltype      = NULL;
131     mtop->nmolblock    = 0;
132     mtop->molblock     = NULL;
133     mtop->maxres_renum = 0;
134     mtop->maxresnr     = -1;
135     init_groups(&mtop->groups);
136     init_block(&mtop->mols);
137     open_symtab(&mtop->symtab);
138 }
139
140 void init_top(t_topology *top)
141 {
142     int i;
143
144     top->name = NULL;
145     init_atom (&(top->atoms));
146     init_atomtypes(&(top->atomtypes));
147     init_block(&top->cgs);
148     init_block(&top->mols);
149     init_blocka(&top->excls);
150     open_symtab(&top->symtab);
151 }
152
153 void init_inputrec(t_inputrec *ir)
154 {
155     memset(ir, 0, (size_t)sizeof(*ir));
156     snew(ir->fepvals, 1);
157     snew(ir->expandedvals, 1);
158     snew(ir->simtempvals, 1);
159 }
160
161 void done_moltype(gmx_moltype_t *molt)
162 {
163     int f;
164
165     done_atom(&molt->atoms);
166     done_block(&molt->cgs);
167     done_blocka(&molt->excls);
168
169     for (f = 0; f < F_NRE; f++)
170     {
171         sfree(molt->ilist[f].iatoms);
172         molt->ilist[f].nalloc = 0;
173     }
174 }
175
176 void done_molblock(gmx_molblock_t *molb)
177 {
178     if (molb->nposres_xA > 0)
179     {
180         molb->nposres_xA = 0;
181         sfree(molb->posres_xA);
182     }
183     if (molb->nposres_xB > 0)
184     {
185         molb->nposres_xB = 0;
186         sfree(molb->posres_xB);
187     }
188 }
189
190 void done_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab)
191 {
192     int i;
193
194     if (bDoneSymtab)
195     {
196         done_symtab(&mtop->symtab);
197     }
198
199     sfree(mtop->ffparams.functype);
200     sfree(mtop->ffparams.iparams);
201
202     for (i = 0; i < mtop->nmoltype; i++)
203     {
204         done_moltype(&mtop->moltype[i]);
205     }
206     sfree(mtop->moltype);
207     for (i = 0; i < mtop->nmolblock; i++)
208     {
209         done_molblock(&mtop->molblock[i]);
210     }
211     sfree(mtop->molblock);
212     done_block(&mtop->mols);
213 }
214
215 void done_top(t_topology *top)
216 {
217     int f;
218
219     sfree(top->idef.functype);
220     sfree(top->idef.iparams);
221     for (f = 0; f < F_NRE; ++f)
222     {
223         sfree(top->idef.il[f].iatoms);
224         top->idef.il[f].iatoms = NULL;
225         top->idef.il[f].nalloc = 0;
226     }
227
228     done_atom (&(top->atoms));
229
230     /* For GB */
231     done_atomtypes(&(top->atomtypes));
232
233     done_symtab(&(top->symtab));
234     done_block(&(top->cgs));
235     done_block(&(top->mols));
236     done_blocka(&(top->excls));
237 }
238
239 static void done_pull_group(t_pull_group *pgrp)
240 {
241     if (pgrp->nat > 0)
242     {
243         sfree(pgrp->ind);
244         sfree(pgrp->ind_loc);
245         sfree(pgrp->weight);
246         sfree(pgrp->weight_loc);
247     }
248 }
249
250 static void done_pull(t_pull *pull)
251 {
252     int i;
253
254     for (i = 0; i < pull->ngroup+1; i++)
255     {
256         done_pull_group(pull->group);
257         done_pull_group(pull->dyna);
258     }
259 }
260
261 void done_inputrec(t_inputrec *ir)
262 {
263     int m;
264
265     for (m = 0; (m < DIM); m++)
266     {
267         if (ir->ex[m].a)
268         {
269             sfree(ir->ex[m].a);
270         }
271         if (ir->ex[m].phi)
272         {
273             sfree(ir->ex[m].phi);
274         }
275         if (ir->et[m].a)
276         {
277             sfree(ir->et[m].a);
278         }
279         if (ir->et[m].phi)
280         {
281             sfree(ir->et[m].phi);
282         }
283     }
284
285     sfree(ir->opts.nrdf);
286     sfree(ir->opts.ref_t);
287     sfree(ir->opts.annealing);
288     sfree(ir->opts.anneal_npoints);
289     sfree(ir->opts.anneal_time);
290     sfree(ir->opts.anneal_temp);
291     sfree(ir->opts.tau_t);
292     sfree(ir->opts.acc);
293     sfree(ir->opts.nFreeze);
294     sfree(ir->opts.QMmethod);
295     sfree(ir->opts.QMbasis);
296     sfree(ir->opts.QMcharge);
297     sfree(ir->opts.QMmult);
298     sfree(ir->opts.bSH);
299     sfree(ir->opts.CASorbitals);
300     sfree(ir->opts.CASelectrons);
301     sfree(ir->opts.SAon);
302     sfree(ir->opts.SAoff);
303     sfree(ir->opts.SAsteps);
304     sfree(ir->opts.bOPT);
305     sfree(ir->opts.bTS);
306
307     if (ir->pull)
308     {
309         done_pull(ir->pull);
310         sfree(ir->pull);
311     }
312 }
313
314 static void zero_history(history_t *hist)
315 {
316     hist->disre_initf  = 0;
317     hist->ndisrepairs  = 0;
318     hist->disre_rm3tav = NULL;
319     hist->orire_initf  = 0;
320     hist->norire_Dtav  = 0;
321     hist->orire_Dtav   = NULL;
322 }
323
324 static void zero_ekinstate(ekinstate_t *eks)
325 {
326     eks->ekin_n         = 0;
327     eks->ekinh          = NULL;
328     eks->ekinf          = NULL;
329     eks->ekinh_old      = NULL;
330     eks->ekinscalef_nhc = NULL;
331     eks->ekinscaleh_nhc = NULL;
332     eks->vscale_nhc     = NULL;
333     eks->dekindl        = 0;
334     eks->mvcos          = 0;
335 }
336
337 static void init_swapstate(swapstate_t *swapstate)
338 {
339     int ii, ic;
340
341     swapstate->eSwapCoords = 0;
342     swapstate->nAverage    = 0;
343
344     /* Ion/water position swapping */
345     for (ic = 0; ic < eCompNR; ic++)
346     {
347         for (ii = 0; ii < eIonNR; ii++)
348         {
349             swapstate->nat_req[ic][ii]        = 0;
350             swapstate->nat_req_p[ic][ii]      = NULL;
351             swapstate->inflow_netto[ic][ii]   = 0;
352             swapstate->inflow_netto_p[ic][ii] = NULL;
353             swapstate->nat_past[ic][ii]       = NULL;
354             swapstate->nat_past_p[ic][ii]     = NULL;
355             swapstate->fluxfromAtoB[ic][ii]   = 0;
356             swapstate->fluxfromAtoB_p[ic][ii] = NULL;
357         }
358     }
359     swapstate->fluxleak               = NULL;
360     swapstate->nions                  = 0;
361     swapstate->comp_from              = NULL;
362     swapstate->channel_label          = NULL;
363     swapstate->bFromCpt               = 0;
364     swapstate->nat[eChan0]            = 0;
365     swapstate->nat[eChan1]            = 0;
366     swapstate->xc_old_whole[eChan0]   = NULL;
367     swapstate->xc_old_whole[eChan1]   = NULL;
368     swapstate->xc_old_whole_p[eChan0] = NULL;
369     swapstate->xc_old_whole_p[eChan1] = NULL;
370 }
371
372 void init_energyhistory(energyhistory_t * enerhist)
373 {
374     enerhist->nener = 0;
375
376     enerhist->ener_ave     = NULL;
377     enerhist->ener_sum     = NULL;
378     enerhist->ener_sum_sim = NULL;
379     enerhist->dht          = NULL;
380
381     enerhist->nsteps     = 0;
382     enerhist->nsum       = 0;
383     enerhist->nsteps_sim = 0;
384     enerhist->nsum_sim   = 0;
385
386     enerhist->dht = NULL;
387 }
388
389 static void done_delta_h_history(delta_h_history_t *dht)
390 {
391     int i;
392
393     for (i = 0; i < dht->nndh; i++)
394     {
395         sfree(dht->dh[i]);
396     }
397     sfree(dht->dh);
398     sfree(dht->ndh);
399 }
400
401 void done_energyhistory(energyhistory_t * enerhist)
402 {
403     sfree(enerhist->ener_ave);
404     sfree(enerhist->ener_sum);
405     sfree(enerhist->ener_sum_sim);
406
407     if (enerhist->dht != NULL)
408     {
409         done_delta_h_history(enerhist->dht);
410         sfree(enerhist->dht);
411     }
412 }
413
414 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
415 {
416     int i, j;
417
418     state->ngtc          = ngtc;
419     state->nnhpres       = nnhpres;
420     state->nhchainlength = nhchainlength;
421     if (state->ngtc > 0)
422     {
423         snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
424         snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
425         snew(state->therm_integral, state->ngtc);
426         for (i = 0; i < state->ngtc; i++)
427         {
428             for (j = 0; j < state->nhchainlength; j++)
429             {
430                 state->nosehoover_xi[i*state->nhchainlength + j]   = 0.0;
431                 state->nosehoover_vxi[i*state->nhchainlength + j]  = 0.0;
432             }
433         }
434         for (i = 0; i < state->ngtc; i++)
435         {
436             state->therm_integral[i]  = 0.0;
437         }
438     }
439     else
440     {
441         state->nosehoover_xi  = NULL;
442         state->nosehoover_vxi = NULL;
443         state->therm_integral = NULL;
444     }
445
446     if (state->nnhpres > 0)
447     {
448         snew(state->nhpres_xi, state->nhchainlength*nnhpres);
449         snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
450         for (i = 0; i < nnhpres; i++)
451         {
452             for (j = 0; j < state->nhchainlength; j++)
453             {
454                 state->nhpres_xi[i*nhchainlength + j]   = 0.0;
455                 state->nhpres_vxi[i*nhchainlength + j]  = 0.0;
456             }
457         }
458     }
459     else
460     {
461         state->nhpres_xi  = NULL;
462         state->nhpres_vxi = NULL;
463     }
464 }
465
466
467 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
468 {
469     int i;
470
471     state->natoms = natoms;
472     state->flags  = 0;
473     state->lambda = 0;
474     snew(state->lambda, efptNR);
475     for (i = 0; i < efptNR; i++)
476     {
477         state->lambda[i] = 0;
478     }
479     state->veta   = 0;
480     clear_mat(state->box);
481     clear_mat(state->box_rel);
482     clear_mat(state->boxv);
483     clear_mat(state->pres_prev);
484     clear_mat(state->svir_prev);
485     clear_mat(state->fvir_prev);
486     init_gtc_state(state, ngtc, nnhpres, nhchainlength);
487     state->nalloc = state->natoms;
488     if (state->nalloc > 0)
489     {
490         snew(state->x, state->nalloc);
491         snew(state->v, state->nalloc);
492     }
493     else
494     {
495         state->x = NULL;
496         state->v = NULL;
497     }
498     state->sd_X = NULL;
499     state->cg_p = NULL;
500     zero_history(&state->hist);
501     zero_ekinstate(&state->ekinstate);
502     init_energyhistory(&state->enerhist);
503     init_df_history(&state->dfhist, nlambda);
504     init_swapstate(&state->swapstate);
505     state->ddp_count       = 0;
506     state->ddp_count_cg_gl = 0;
507     state->cg_gl           = NULL;
508     state->cg_gl_nalloc    = 0;
509 }
510
511 void done_state(t_state *state)
512 {
513     if (state->x)
514     {
515         sfree(state->x);
516     }
517     if (state->v)
518     {
519         sfree(state->v);
520     }
521     if (state->sd_X)
522     {
523         sfree(state->sd_X);
524     }
525     if (state->cg_p)
526     {
527         sfree(state->cg_p);
528     }
529     state->nalloc = 0;
530     if (state->cg_gl)
531     {
532         sfree(state->cg_gl);
533     }
534     state->cg_gl_nalloc = 0;
535     if (state->lambda)
536     {
537         sfree(state->lambda);
538     }
539     if (state->ngtc > 0)
540     {
541         sfree(state->nosehoover_xi);
542         sfree(state->nosehoover_vxi);
543         sfree(state->therm_integral);
544     }
545 }
546
547 t_state *serial_init_local_state(t_state *state_global)
548 {
549     int      i;
550     t_state *state_local;
551
552     snew(state_local, 1);
553
554     /* Copy all the contents */
555     *state_local = *state_global;
556     snew(state_local->lambda, efptNR);
557     /* local storage for lambda */
558     for (i = 0; i < efptNR; i++)
559     {
560         state_local->lambda[i] = state_global->lambda[i];
561     }
562
563     return state_local;
564 }
565
566 static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
567 {
568     int d, d2;
569
570     for (d = YY; d <= ZZ; d++)
571     {
572         for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
573         {
574             /* We need to check if this box component is deformed
575              * or if deformation of another component might cause
576              * changes in this component due to box corrections.
577              */
578             if (ir->deform[d][d2] == 0 &&
579                 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
580                   (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
581             {
582                 if (bInit)
583                 {
584                     box_rel[d][d2] = b[d][d2]/b[XX][XX];
585                 }
586                 else
587                 {
588                     b[d][d2] = b[XX][XX]*box_rel[d][d2];
589                 }
590             }
591         }
592     }
593 }
594
595 void set_box_rel(t_inputrec *ir, t_state *state)
596 {
597     /* Make sure the box obeys the restrictions before we fix the ratios */
598     correct_box(NULL, 0, state->box, NULL);
599
600     clear_mat(state->box_rel);
601
602     if (PRESERVE_SHAPE(*ir))
603     {
604         do_box_rel(ir, state->box_rel, state->box, TRUE);
605     }
606 }
607
608 void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
609 {
610     if (PRESERVE_SHAPE(*ir))
611     {
612         do_box_rel(ir, box_rel, b, FALSE);
613     }
614 }
615
616 real max_cutoff(real cutoff1, real cutoff2)
617 {
618     if (cutoff1 == 0 || cutoff2 == 0)
619     {
620         return 0;
621     }
622     else
623     {
624         return max(cutoff1, cutoff2);
625     }
626 }
627
628 void init_df_history(df_history_t *dfhist, int nlambda)
629 {
630     int i;
631
632     dfhist->nlambda  = nlambda;
633     dfhist->bEquil   = 0;
634     dfhist->wl_delta = 0;
635
636     if (nlambda > 0)
637     {
638         snew(dfhist->sum_weights, dfhist->nlambda);
639         snew(dfhist->sum_dg, dfhist->nlambda);
640         snew(dfhist->sum_minvar, dfhist->nlambda);
641         snew(dfhist->sum_variance, dfhist->nlambda);
642         snew(dfhist->n_at_lam, dfhist->nlambda);
643         snew(dfhist->wl_histo, dfhist->nlambda);
644
645         /* allocate transition matrices here */
646         snew(dfhist->Tij, dfhist->nlambda);
647         snew(dfhist->Tij_empirical, dfhist->nlambda);
648
649         /* allocate accumulators for various transition matrix
650            free energy methods here */
651         snew(dfhist->accum_p, dfhist->nlambda);
652         snew(dfhist->accum_m, dfhist->nlambda);
653         snew(dfhist->accum_p2, dfhist->nlambda);
654         snew(dfhist->accum_m2, dfhist->nlambda);
655
656         for (i = 0; i < dfhist->nlambda; i++)
657         {
658             snew(dfhist->Tij[i], dfhist->nlambda);
659             snew(dfhist->Tij_empirical[i], dfhist->nlambda);
660             snew((dfhist->accum_p)[i], dfhist->nlambda);
661             snew((dfhist->accum_m)[i], dfhist->nlambda);
662             snew((dfhist->accum_p2)[i], dfhist->nlambda);
663             snew((dfhist->accum_m2)[i], dfhist->nlambda);
664         }
665     }
666 }
667
668 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
669 {
670     int i, j;
671
672     /* Currently, there should not be any difference in nlambda between the two,
673        but this is included for completeness for potential later functionality */
674     df_dest->nlambda  = df_source->nlambda;
675     df_dest->bEquil   = df_source->bEquil;
676     df_dest->wl_delta = df_source->wl_delta;
677
678     for (i = 0; i < df_dest->nlambda; i++)
679     {
680         df_dest->sum_weights[i]  = df_source->sum_weights[i];
681         df_dest->sum_dg[i]       = df_source->sum_dg[i];
682         df_dest->sum_minvar[i]   = df_source->sum_minvar[i];
683         df_dest->sum_variance[i] = df_source->sum_variance[i];
684         df_dest->n_at_lam[i]     = df_source->n_at_lam[i];
685         df_dest->wl_histo[i]     = df_source->wl_histo[i];
686     }
687
688     for (i = 0; i < df_dest->nlambda; i++)
689     {
690         for (j = 0; j < df_dest->nlambda; j++)
691         {
692             df_dest->accum_p[i][j]        = df_source->accum_p[i][j];
693             df_dest->accum_m[i][j]        = df_source->accum_m[i][j];
694             df_dest->accum_p2[i][j]       = df_source->accum_p2[i][j];
695             df_dest->accum_m2[i][j]       = df_source->accum_m2[i][j];
696             df_dest->Tij[i][j]            = df_source->Tij[i][j];
697             df_dest->Tij_empirical[i][j]  = df_source->Tij_empirical[i][j];
698         }
699     }
700 }
701
702 void done_df_history(df_history_t *dfhist)
703 {
704     int i;
705
706     if (dfhist->nlambda > 0)
707     {
708         sfree(dfhist->n_at_lam);
709         sfree(dfhist->wl_histo);
710         sfree(dfhist->sum_weights);
711         sfree(dfhist->sum_dg);
712         sfree(dfhist->sum_minvar);
713         sfree(dfhist->sum_variance);
714
715         for (i = 0; i < dfhist->nlambda; i++)
716         {
717             sfree(dfhist->Tij[i]);
718             sfree(dfhist->Tij_empirical[i]);
719             sfree(dfhist->accum_p[i]);
720             sfree(dfhist->accum_m[i]);
721             sfree(dfhist->accum_p2[i]);
722             sfree(dfhist->accum_m2[i]);
723         }
724     }
725     dfhist->bEquil   = 0;
726     dfhist->nlambda  = 0;
727     dfhist->wl_delta = 0;
728 }