Merge branch 'master' into pygromacs
[alexxy/gromacs.git] / src / gromacs / gmxlib / typedefs.cpp
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,2015, 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 "gmxpre.h"
39
40 #include "gromacs/legacyheaders/typedefs.h"
41
42 #include <cstring>
43
44 #include <algorithm>
45
46 #include "gromacs/legacyheaders/macros.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/pbcutil/pbc.h"
49 #include "gromacs/random/random.h"
50 #include "gromacs/utility/smalloc.h"
51
52 /* The source code in this file should be thread-safe.
53       Please keep it that way. */
54
55 int gmx_int64_to_int(gmx_int64_t step, const char *warn)
56 {
57     int i;
58
59     i = static_cast<int>(step);
60
61     if (warn != NULL && (step < INT_MIN || step > INT_MAX))
62     {
63         fprintf(stderr, "\nWARNING during %s:\n", warn);
64         fprintf(stderr, "step value ");
65         fprintf(stderr, "%" GMX_PRId64, step);
66         fprintf(stderr, " does not fit in int, converted to %d\n\n", i);
67     }
68
69     return i;
70 }
71
72 void init_inputrec(t_inputrec *ir)
73 {
74     std::memset(ir, 0, sizeof(*ir));
75     snew(ir->fepvals, 1);
76     snew(ir->expandedvals, 1);
77     snew(ir->simtempvals, 1);
78 }
79
80 static void done_pull_group(t_pull_group *pgrp)
81 {
82     if (pgrp->nat > 0)
83     {
84         sfree(pgrp->ind);
85         sfree(pgrp->weight);
86     }
87 }
88
89 static void done_pull_params(pull_params_t *pull)
90 {
91     int i;
92
93     for (i = 0; i < pull->ngroup+1; i++)
94     {
95         done_pull_group(pull->group);
96     }
97
98     sfree(pull->group);
99     sfree(pull->coord);
100 }
101
102 void done_inputrec(t_inputrec *ir)
103 {
104     int m;
105
106     for (m = 0; (m < DIM); m++)
107     {
108         if (ir->ex[m].a)
109         {
110             sfree(ir->ex[m].a);
111         }
112         if (ir->ex[m].phi)
113         {
114             sfree(ir->ex[m].phi);
115         }
116         if (ir->et[m].a)
117         {
118             sfree(ir->et[m].a);
119         }
120         if (ir->et[m].phi)
121         {
122             sfree(ir->et[m].phi);
123         }
124     }
125
126     sfree(ir->opts.nrdf);
127     sfree(ir->opts.ref_t);
128     sfree(ir->opts.annealing);
129     sfree(ir->opts.anneal_npoints);
130     sfree(ir->opts.anneal_time);
131     sfree(ir->opts.anneal_temp);
132     sfree(ir->opts.tau_t);
133     sfree(ir->opts.acc);
134     sfree(ir->opts.nFreeze);
135     sfree(ir->opts.QMmethod);
136     sfree(ir->opts.QMbasis);
137     sfree(ir->opts.QMcharge);
138     sfree(ir->opts.QMmult);
139     sfree(ir->opts.bSH);
140     sfree(ir->opts.CASorbitals);
141     sfree(ir->opts.CASelectrons);
142     sfree(ir->opts.SAon);
143     sfree(ir->opts.SAoff);
144     sfree(ir->opts.SAsteps);
145     sfree(ir->opts.bOPT);
146     sfree(ir->opts.bTS);
147
148     if (ir->pull)
149     {
150         done_pull_params(ir->pull);
151         sfree(ir->pull);
152     }
153 }
154
155 static void zero_history(history_t *hist)
156 {
157     hist->disre_initf  = 0;
158     hist->ndisrepairs  = 0;
159     hist->disre_rm3tav = NULL;
160     hist->orire_initf  = 0;
161     hist->norire_Dtav  = 0;
162     hist->orire_Dtav   = NULL;
163 }
164
165 static void zero_ekinstate(ekinstate_t *eks)
166 {
167     eks->ekin_n         = 0;
168     eks->ekinh          = NULL;
169     eks->ekinf          = NULL;
170     eks->ekinh_old      = NULL;
171     eks->ekinscalef_nhc = NULL;
172     eks->ekinscaleh_nhc = NULL;
173     eks->vscale_nhc     = NULL;
174     eks->dekindl        = 0;
175     eks->mvcos          = 0;
176 }
177
178 static void init_swapstate(swapstate_t *swapstate)
179 {
180     int ii, ic;
181
182     swapstate->eSwapCoords = 0;
183     swapstate->nAverage    = 0;
184
185     /* Ion/water position swapping */
186     for (ic = 0; ic < eCompNR; ic++)
187     {
188         for (ii = 0; ii < eIonNR; ii++)
189         {
190             swapstate->nat_req[ic][ii]        = 0;
191             swapstate->nat_req_p[ic][ii]      = NULL;
192             swapstate->inflow_netto[ic][ii]   = 0;
193             swapstate->inflow_netto_p[ic][ii] = NULL;
194             swapstate->nat_past[ic][ii]       = NULL;
195             swapstate->nat_past_p[ic][ii]     = NULL;
196             swapstate->fluxfromAtoB[ic][ii]   = 0;
197             swapstate->fluxfromAtoB_p[ic][ii] = NULL;
198         }
199     }
200     swapstate->fluxleak               = NULL;
201     swapstate->nions                  = 0;
202     swapstate->comp_from              = NULL;
203     swapstate->channel_label          = NULL;
204     swapstate->bFromCpt               = 0;
205     swapstate->nat[eChan0]            = 0;
206     swapstate->nat[eChan1]            = 0;
207     swapstate->xc_old_whole[eChan0]   = NULL;
208     swapstate->xc_old_whole[eChan1]   = NULL;
209     swapstate->xc_old_whole_p[eChan0] = NULL;
210     swapstate->xc_old_whole_p[eChan1] = NULL;
211 }
212
213 void init_energyhistory(energyhistory_t * enerhist)
214 {
215     enerhist->nener = 0;
216
217     enerhist->ener_ave     = NULL;
218     enerhist->ener_sum     = NULL;
219     enerhist->ener_sum_sim = NULL;
220     enerhist->dht          = NULL;
221
222     enerhist->nsteps     = 0;
223     enerhist->nsum       = 0;
224     enerhist->nsteps_sim = 0;
225     enerhist->nsum_sim   = 0;
226
227     enerhist->dht = NULL;
228 }
229
230 static void done_delta_h_history(delta_h_history_t *dht)
231 {
232     int i;
233
234     for (i = 0; i < dht->nndh; i++)
235     {
236         sfree(dht->dh[i]);
237     }
238     sfree(dht->dh);
239     sfree(dht->ndh);
240 }
241
242 void done_energyhistory(energyhistory_t * enerhist)
243 {
244     sfree(enerhist->ener_ave);
245     sfree(enerhist->ener_sum);
246     sfree(enerhist->ener_sum_sim);
247
248     if (enerhist->dht != NULL)
249     {
250         done_delta_h_history(enerhist->dht);
251         sfree(enerhist->dht);
252     }
253 }
254
255 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
256 {
257     int i, j;
258
259     state->ngtc          = ngtc;
260     state->nnhpres       = nnhpres;
261     state->nhchainlength = nhchainlength;
262     if (state->ngtc > 0)
263     {
264         snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
265         snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
266         snew(state->therm_integral, state->ngtc);
267         for (i = 0; i < state->ngtc; i++)
268         {
269             for (j = 0; j < state->nhchainlength; j++)
270             {
271                 state->nosehoover_xi[i*state->nhchainlength + j]   = 0.0;
272                 state->nosehoover_vxi[i*state->nhchainlength + j]  = 0.0;
273             }
274         }
275         for (i = 0; i < state->ngtc; i++)
276         {
277             state->therm_integral[i]  = 0.0;
278         }
279     }
280     else
281     {
282         state->nosehoover_xi  = NULL;
283         state->nosehoover_vxi = NULL;
284         state->therm_integral = NULL;
285     }
286
287     if (state->nnhpres > 0)
288     {
289         snew(state->nhpres_xi, state->nhchainlength*nnhpres);
290         snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
291         for (i = 0; i < nnhpres; i++)
292         {
293             for (j = 0; j < state->nhchainlength; j++)
294             {
295                 state->nhpres_xi[i*nhchainlength + j]   = 0.0;
296                 state->nhpres_vxi[i*nhchainlength + j]  = 0.0;
297             }
298         }
299     }
300     else
301     {
302         state->nhpres_xi  = NULL;
303         state->nhpres_vxi = NULL;
304     }
305 }
306
307
308 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
309 {
310     int i;
311
312     state->natoms    = natoms;
313     state->flags     = 0;
314     state->fep_state = 0;
315     state->lambda    = 0;
316     snew(state->lambda, efptNR);
317     for (i = 0; i < efptNR; i++)
318     {
319         state->lambda[i] = 0;
320     }
321     state->veta   = 0;
322     clear_mat(state->box);
323     clear_mat(state->box_rel);
324     clear_mat(state->boxv);
325     clear_mat(state->pres_prev);
326     clear_mat(state->svir_prev);
327     clear_mat(state->fvir_prev);
328     init_gtc_state(state, ngtc, nnhpres, nhchainlength);
329     state->nalloc = state->natoms;
330     if (state->nalloc > 0)
331     {
332         snew(state->x, state->nalloc);
333         snew(state->v, state->nalloc);
334     }
335     else
336     {
337         state->x = NULL;
338         state->v = NULL;
339     }
340     state->sd_X = NULL;
341     state->cg_p = NULL;
342     zero_history(&state->hist);
343     zero_ekinstate(&state->ekinstate);
344     init_energyhistory(&state->enerhist);
345     init_df_history(&state->dfhist, nlambda);
346     init_swapstate(&state->swapstate);
347     state->ddp_count       = 0;
348     state->ddp_count_cg_gl = 0;
349     state->cg_gl           = NULL;
350     state->cg_gl_nalloc    = 0;
351 }
352
353 void done_state(t_state *state)
354 {
355     if (state->x)
356     {
357         sfree(state->x);
358     }
359     if (state->v)
360     {
361         sfree(state->v);
362     }
363     if (state->sd_X)
364     {
365         sfree(state->sd_X);
366     }
367     if (state->cg_p)
368     {
369         sfree(state->cg_p);
370     }
371     state->nalloc = 0;
372     if (state->cg_gl)
373     {
374         sfree(state->cg_gl);
375     }
376     state->cg_gl_nalloc = 0;
377     if (state->lambda)
378     {
379         sfree(state->lambda);
380     }
381     if (state->ngtc > 0)
382     {
383         sfree(state->nosehoover_xi);
384         sfree(state->nosehoover_vxi);
385         sfree(state->therm_integral);
386     }
387 }
388
389 t_state *serial_init_local_state(t_state *state_global)
390 {
391     int      i;
392     t_state *state_local;
393
394     snew(state_local, 1);
395
396     /* Copy all the contents */
397     *state_local = *state_global;
398     snew(state_local->lambda, efptNR);
399     /* local storage for lambda */
400     for (i = 0; i < efptNR; i++)
401     {
402         state_local->lambda[i] = state_global->lambda[i];
403     }
404
405     return state_local;
406 }
407
408 static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
409 {
410     int d, d2;
411
412     for (d = YY; d <= ZZ; d++)
413     {
414         for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
415         {
416             /* We need to check if this box component is deformed
417              * or if deformation of another component might cause
418              * changes in this component due to box corrections.
419              */
420             if (ir->deform[d][d2] == 0 &&
421                 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
422                   (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
423             {
424                 if (bInit)
425                 {
426                     box_rel[d][d2] = b[d][d2]/b[XX][XX];
427                 }
428                 else
429                 {
430                     b[d][d2] = b[XX][XX]*box_rel[d][d2];
431                 }
432             }
433         }
434     }
435 }
436
437 void set_box_rel(t_inputrec *ir, t_state *state)
438 {
439     /* Make sure the box obeys the restrictions before we fix the ratios */
440     correct_box(NULL, 0, state->box, NULL);
441
442     clear_mat(state->box_rel);
443
444     if (PRESERVE_SHAPE(*ir))
445     {
446         do_box_rel(ir, state->box_rel, state->box, TRUE);
447     }
448 }
449
450 void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
451 {
452     if (PRESERVE_SHAPE(*ir))
453     {
454         do_box_rel(ir, box_rel, b, FALSE);
455     }
456 }
457
458 real max_cutoff(real cutoff1, real cutoff2)
459 {
460     if (cutoff1 == 0 || cutoff2 == 0)
461     {
462         return 0;
463     }
464     else
465     {
466         return std::max(cutoff1, cutoff2);
467     }
468 }
469
470 void init_df_history(df_history_t *dfhist, int nlambda)
471 {
472     int i;
473
474     dfhist->nlambda  = nlambda;
475     dfhist->bEquil   = 0;
476     dfhist->wl_delta = 0;
477
478     if (nlambda > 0)
479     {
480         snew(dfhist->sum_weights, dfhist->nlambda);
481         snew(dfhist->sum_dg, dfhist->nlambda);
482         snew(dfhist->sum_minvar, dfhist->nlambda);
483         snew(dfhist->sum_variance, dfhist->nlambda);
484         snew(dfhist->n_at_lam, dfhist->nlambda);
485         snew(dfhist->wl_histo, dfhist->nlambda);
486
487         /* allocate transition matrices here */
488         snew(dfhist->Tij, dfhist->nlambda);
489         snew(dfhist->Tij_empirical, dfhist->nlambda);
490
491         /* allocate accumulators for various transition matrix
492            free energy methods here */
493         snew(dfhist->accum_p, dfhist->nlambda);
494         snew(dfhist->accum_m, dfhist->nlambda);
495         snew(dfhist->accum_p2, dfhist->nlambda);
496         snew(dfhist->accum_m2, dfhist->nlambda);
497
498         for (i = 0; i < dfhist->nlambda; i++)
499         {
500             snew(dfhist->Tij[i], dfhist->nlambda);
501             snew(dfhist->Tij_empirical[i], dfhist->nlambda);
502             snew((dfhist->accum_p)[i], dfhist->nlambda);
503             snew((dfhist->accum_m)[i], dfhist->nlambda);
504             snew((dfhist->accum_p2)[i], dfhist->nlambda);
505             snew((dfhist->accum_m2)[i], dfhist->nlambda);
506         }
507     }
508 }
509
510 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
511 {
512     int i, j;
513
514     /* Currently, there should not be any difference in nlambda between the two,
515        but this is included for completeness for potential later functionality */
516     df_dest->nlambda  = df_source->nlambda;
517     df_dest->bEquil   = df_source->bEquil;
518     df_dest->wl_delta = df_source->wl_delta;
519
520     for (i = 0; i < df_dest->nlambda; i++)
521     {
522         df_dest->sum_weights[i]  = df_source->sum_weights[i];
523         df_dest->sum_dg[i]       = df_source->sum_dg[i];
524         df_dest->sum_minvar[i]   = df_source->sum_minvar[i];
525         df_dest->sum_variance[i] = df_source->sum_variance[i];
526         df_dest->n_at_lam[i]     = df_source->n_at_lam[i];
527         df_dest->wl_histo[i]     = df_source->wl_histo[i];
528     }
529
530     for (i = 0; i < df_dest->nlambda; i++)
531     {
532         for (j = 0; j < df_dest->nlambda; j++)
533         {
534             df_dest->accum_p[i][j]        = df_source->accum_p[i][j];
535             df_dest->accum_m[i][j]        = df_source->accum_m[i][j];
536             df_dest->accum_p2[i][j]       = df_source->accum_p2[i][j];
537             df_dest->accum_m2[i][j]       = df_source->accum_m2[i][j];
538             df_dest->Tij[i][j]            = df_source->Tij[i][j];
539             df_dest->Tij_empirical[i][j]  = df_source->Tij_empirical[i][j];
540         }
541     }
542 }
543
544 void done_df_history(df_history_t *dfhist)
545 {
546     int i;
547
548     if (dfhist->nlambda > 0)
549     {
550         sfree(dfhist->n_at_lam);
551         sfree(dfhist->wl_histo);
552         sfree(dfhist->sum_weights);
553         sfree(dfhist->sum_dg);
554         sfree(dfhist->sum_minvar);
555         sfree(dfhist->sum_variance);
556
557         for (i = 0; i < dfhist->nlambda; i++)
558         {
559             sfree(dfhist->Tij[i]);
560             sfree(dfhist->Tij_empirical[i]);
561             sfree(dfhist->accum_p[i]);
562             sfree(dfhist->accum_m[i]);
563             sfree(dfhist->accum_p2[i]);
564             sfree(dfhist->accum_m2[i]);
565         }
566     }
567     dfhist->bEquil   = 0;
568     dfhist->nlambda  = 0;
569     dfhist->wl_delta = 0;
570 }