1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROningen Mixture of Alchemy and Childrens' Stories
36 /* This file is completely threadsafe - keep it that way! */
48 #include "thread_mpi.h"
51 /* The source code in this file should be thread-safe.
52 Please keep it that way. */
56 static gmx_bool bOverAllocDD
=FALSE
;
58 static tMPI_Thread_mutex_t over_alloc_mutex
=TMPI_THREAD_MUTEX_INITIALIZER
;
62 void set_over_alloc_dd(gmx_bool set
)
65 tMPI_Thread_mutex_lock(&over_alloc_mutex
);
66 /* we just make sure that we don't set this at the same time.
67 We don't worry too much about reading this rarely-set variable */
71 tMPI_Thread_mutex_unlock(&over_alloc_mutex
);
75 int over_alloc_dd(int n
)
78 return OVER_ALLOC_FAC
*n
+ 100;
83 int gmx_large_int_to_int(gmx_large_int_t step
,const char *warn
)
89 if (warn
!= NULL
&& (step
< INT_MIN
|| step
> INT_MAX
)) {
90 fprintf(stderr
,"\nWARNING during %s:\n",warn
);
91 fprintf(stderr
,"step value ");
92 fprintf(stderr
,gmx_large_int_pfmt
,step
);
93 fprintf(stderr
," does not fit in int, converted to %d\n\n",i
);
99 char *gmx_step_str(gmx_large_int_t i
,char *buf
)
101 sprintf(buf
,gmx_large_int_pfmt
,i
);
106 void init_block(t_block
*block
)
111 block
->nalloc_index
= 1;
112 snew(block
->index
,block
->nalloc_index
);
116 void init_blocka(t_blocka
*block
)
122 block
->nalloc_index
= 1;
123 snew(block
->index
,block
->nalloc_index
);
129 void init_atom(t_atoms
*at
)
143 void init_atomtypes(t_atomtypes
*at
)
148 at
->atomnumber
= NULL
;
149 at
->gb_radius
= NULL
;
153 void init_groups(gmx_groups_t
*groups
)
157 groups
->ngrpname
= 0;
158 groups
->grpname
= NULL
;
159 for(g
=0; (g
<egcNR
); g
++) {
160 groups
->grps
[g
].nm_ind
= NULL
;
161 groups
->ngrpnr
[g
] = 0;
162 groups
->grpnr
[g
] = NULL
;
167 void init_mtop(gmx_mtop_t
*mtop
)
171 mtop
->moltype
= NULL
;
173 mtop
->molblock
= NULL
;
174 mtop
->maxres_renum
= 0;
176 init_groups(&mtop
->groups
);
177 init_block(&mtop
->mols
);
178 open_symtab(&mtop
->symtab
);
181 void init_top (t_topology
*top
)
186 init_atom (&(top
->atoms
));
187 init_atomtypes(&(top
->atomtypes
));
188 init_block(&top
->cgs
);
189 init_block(&top
->mols
);
190 init_blocka(&top
->excls
);
191 open_symtab(&top
->symtab
);
194 void init_inputrec(t_inputrec
*ir
)
196 memset(ir
,0,(size_t)sizeof(*ir
));
199 void stupid_fill_block(t_block
*grp
,int natom
,gmx_bool bOneIndexGroup
)
203 if (bOneIndexGroup
) {
204 grp
->nalloc_index
= 2;
205 snew(grp
->index
,grp
->nalloc_index
);
211 grp
->nalloc_index
= natom
+1;
212 snew(grp
->index
,grp
->nalloc_index
);
213 snew(grp
->index
,natom
+1);
214 for(i
=0; (i
<=natom
); i
++)
220 void stupid_fill_blocka(t_blocka
*grp
,int natom
)
224 grp
->nalloc_a
= natom
;
225 snew(grp
->a
,grp
->nalloc_a
);
226 for(i
=0; (i
<natom
); i
++)
230 grp
->nalloc_index
= natom
+ 1;
231 snew(grp
->index
,grp
->nalloc_index
);
232 for(i
=0; (i
<=natom
); i
++)
237 void copy_blocka(const t_blocka
*src
,t_blocka
*dest
)
242 dest
->nalloc_index
= dest
->nr
+ 1;
243 snew(dest
->index
,dest
->nalloc_index
);
244 for(i
=0; i
<dest
->nr
+1; i
++) {
245 dest
->index
[i
] = src
->index
[i
];
247 dest
->nra
= src
->nra
;
248 dest
->nalloc_a
= dest
->nra
+ 1;
249 snew(dest
->a
,dest
->nalloc_a
);
250 for(i
=0; i
<dest
->nra
+1; i
++) {
251 dest
->a
[i
] = src
->a
[i
];
255 void done_block(t_block
*block
)
259 block
->nalloc_index
= 0;
262 void done_blocka(t_blocka
*block
)
269 block
->nalloc_index
= 0;
273 void done_atom (t_atoms
*at
)
281 sfree(at
->atomtypeB
);
284 void done_atomtypes(t_atomtypes
*atype
)
287 sfree(atype
->radius
);
289 sfree(atype
->surftens
);
290 sfree(atype
->atomnumber
);
291 sfree(atype
->gb_radius
);
295 void done_moltype(gmx_moltype_t
*molt
)
299 done_atom(&molt
->atoms
);
300 done_block(&molt
->cgs
);
301 done_blocka(&molt
->excls
);
303 for(f
=0; f
<F_NRE
; f
++) {
304 sfree(molt
->ilist
[f
].iatoms
);
305 molt
->ilist
[f
].nalloc
= 0;
309 void done_molblock(gmx_molblock_t
*molb
)
311 if (molb
->nposres_xA
> 0) {
312 molb
->nposres_xA
= 0;
313 free(molb
->posres_xA
);
315 if (molb
->nposres_xB
> 0) {
316 molb
->nposres_xB
= 0;
317 free(molb
->posres_xB
);
321 void done_mtop(gmx_mtop_t
*mtop
,gmx_bool bDoneSymtab
)
326 done_symtab(&mtop
->symtab
);
329 sfree(mtop
->ffparams
.functype
);
330 sfree(mtop
->ffparams
.iparams
);
332 for(i
=0; i
<mtop
->nmoltype
; i
++) {
333 done_moltype(&mtop
->moltype
[i
]);
335 sfree(mtop
->moltype
);
336 for(i
=0; i
<mtop
->nmolblock
; i
++) {
337 done_molblock(&mtop
->molblock
[i
]);
339 sfree(mtop
->molblock
);
340 done_block(&mtop
->mols
);
343 void done_top(t_topology
*top
)
347 sfree(top
->idef
.functype
);
348 sfree(top
->idef
.iparams
);
349 for (f
= 0; f
< F_NRE
; ++f
)
351 sfree(top
->idef
.il
[f
].iatoms
);
352 top
->idef
.il
[f
].iatoms
= NULL
;
353 top
->idef
.il
[f
].nalloc
= 0;
356 done_atom (&(top
->atoms
));
359 done_atomtypes(&(top
->atomtypes
));
361 done_symtab(&(top
->symtab
));
362 done_block(&(top
->cgs
));
363 done_block(&(top
->mols
));
364 done_blocka(&(top
->excls
));
367 static void done_pullgrp(t_pullgrp
*pgrp
)
370 sfree(pgrp
->ind_loc
);
372 sfree(pgrp
->weight_loc
);
375 static void done_pull(t_pull
*pull
)
379 for(i
=0; i
<pull
->ngrp
+1; i
++) {
380 done_pullgrp(pull
->grp
);
381 done_pullgrp(pull
->dyna
);
385 void done_inputrec(t_inputrec
*ir
)
389 for(m
=0; (m
<DIM
); m
++) {
390 if (ir
->ex
[m
].a
) sfree(ir
->ex
[m
].a
);
391 if (ir
->ex
[m
].phi
) sfree(ir
->ex
[m
].phi
);
392 if (ir
->et
[m
].a
) sfree(ir
->et
[m
].a
);
393 if (ir
->et
[m
].phi
) sfree(ir
->et
[m
].phi
);
396 sfree(ir
->opts
.nrdf
);
397 sfree(ir
->opts
.ref_t
);
398 sfree(ir
->opts
.annealing
);
399 sfree(ir
->opts
.anneal_npoints
);
400 sfree(ir
->opts
.anneal_time
);
401 sfree(ir
->opts
.anneal_temp
);
402 sfree(ir
->opts
.tau_t
);
404 sfree(ir
->opts
.nFreeze
);
405 sfree(ir
->opts
.QMmethod
);
406 sfree(ir
->opts
.QMbasis
);
407 sfree(ir
->opts
.QMcharge
);
408 sfree(ir
->opts
.QMmult
);
410 sfree(ir
->opts
.CASorbitals
);
411 sfree(ir
->opts
.CASelectrons
);
412 sfree(ir
->opts
.SAon
);
413 sfree(ir
->opts
.SAoff
);
414 sfree(ir
->opts
.SAsteps
);
415 sfree(ir
->opts
.bOPT
);
424 static void zero_ekinstate(ekinstate_t
*eks
)
429 eks
->ekinh_old
= NULL
;
430 eks
->ekinscalef_nhc
= NULL
;
431 eks
->ekinscaleh_nhc
= NULL
;
432 eks
->vscale_nhc
= NULL
;
437 void init_energyhistory(energyhistory_t
* enerhist
)
441 enerhist
->ener_ave
= NULL
;
442 enerhist
->ener_sum
= NULL
;
443 enerhist
->ener_sum_sim
= NULL
;
444 enerhist
->dht
= NULL
;
446 enerhist
->nsteps
= 0;
448 enerhist
->nsteps_sim
= 0;
449 enerhist
->nsum_sim
= 0;
451 enerhist
->dht
= NULL
;
454 static void done_delta_h_history(delta_h_history_t
*dht
)
458 for(i
=0; i
<dht
->nndh
; i
++)
466 void done_energyhistory(energyhistory_t
* enerhist
)
468 sfree(enerhist
->ener_ave
);
469 sfree(enerhist
->ener_sum
);
470 sfree(enerhist
->ener_sum_sim
);
472 if (enerhist
->dht
!= NULL
)
474 done_delta_h_history(enerhist
->dht
);
475 sfree(enerhist
->dht
);
479 void init_gtc_state(t_state
*state
, int ngtc
, int nnhpres
, int nhchainlength
)
484 state
->nnhpres
= nnhpres
;
485 state
->nhchainlength
= nhchainlength
;
488 snew(state
->nosehoover_xi
,state
->nhchainlength
*state
->ngtc
);
489 snew(state
->nosehoover_vxi
,state
->nhchainlength
*state
->ngtc
);
490 snew(state
->therm_integral
,state
->ngtc
);
491 for(i
=0; i
<state
->ngtc
; i
++)
493 for (j
=0;j
<state
->nhchainlength
;j
++)
495 state
->nosehoover_xi
[i
*state
->nhchainlength
+ j
] = 0.0;
496 state
->nosehoover_vxi
[i
*state
->nhchainlength
+ j
] = 0.0;
499 for(i
=0; i
<state
->ngtc
; i
++) {
500 state
->therm_integral
[i
] = 0.0;
505 state
->nosehoover_xi
= NULL
;
506 state
->nosehoover_vxi
= NULL
;
507 state
->therm_integral
= NULL
;
510 if (state
->nnhpres
> 0)
512 snew(state
->nhpres_xi
,state
->nhchainlength
*nnhpres
);
513 snew(state
->nhpres_vxi
,state
->nhchainlength
*nnhpres
);
514 for(i
=0; i
<nnhpres
; i
++)
516 for (j
=0;j
<state
->nhchainlength
;j
++)
518 state
->nhpres_xi
[i
*nhchainlength
+ j
] = 0.0;
519 state
->nhpres_vxi
[i
*nhchainlength
+ j
] = 0.0;
525 state
->nhpres_xi
= NULL
;
526 state
->nhpres_vxi
= NULL
;
531 void init_state(t_state
*state
, int natoms
, int ngtc
, int nnhpres
, int nhchainlength
)
535 state
->natoms
= natoms
;
540 clear_mat(state
->box
);
541 clear_mat(state
->box_rel
);
542 clear_mat(state
->boxv
);
543 clear_mat(state
->pres_prev
);
544 clear_mat(state
->svir_prev
);
545 clear_mat(state
->fvir_prev
);
546 init_gtc_state(state
,ngtc
,nnhpres
,nhchainlength
);
547 state
->nalloc
= state
->natoms
;
548 if (state
->nalloc
> 0) {
549 snew(state
->x
,state
->nalloc
);
550 snew(state
->v
,state
->nalloc
);
558 zero_ekinstate(&state
->ekinstate
);
560 init_energyhistory(&state
->enerhist
);
562 state
->ddp_count
= 0;
563 state
->ddp_count_cg_gl
= 0;
565 state
->cg_gl_nalloc
= 0;
568 void done_state(t_state
*state
)
570 if (state
->nosehoover_xi
) sfree(state
->nosehoover_xi
);
571 if (state
->x
) sfree(state
->x
);
572 if (state
->v
) sfree(state
->v
);
573 if (state
->sd_X
) sfree(state
->sd_X
);
574 if (state
->cg_p
) sfree(state
->cg_p
);
576 if (state
->cg_gl
) sfree(state
->cg_gl
);
577 state
->cg_gl_nalloc
= 0;
580 static void do_box_rel(t_inputrec
*ir
,matrix box_rel
,matrix b
,gmx_bool bInit
)
584 for(d
=YY
; d
<=ZZ
; d
++) {
585 for(d2
=XX
; d2
<=(ir
->epct
==epctSEMIISOTROPIC
? YY
: ZZ
); d2
++) {
586 /* We need to check if this box component is deformed
587 * or if deformation of another component might cause
588 * changes in this component due to box corrections.
590 if (ir
->deform
[d
][d2
] == 0 &&
591 !(d
== ZZ
&& d2
== XX
&& ir
->deform
[d
][YY
] != 0 &&
592 (b
[YY
][d2
] != 0 || ir
->deform
[YY
][d2
] != 0))) {
594 box_rel
[d
][d2
] = b
[d
][d2
]/b
[XX
][XX
];
596 b
[d
][d2
] = b
[XX
][XX
]*box_rel
[d
][d2
];
603 void set_box_rel(t_inputrec
*ir
,t_state
*state
)
605 /* Make sure the box obeys the restrictions before we fix the ratios */
606 correct_box(NULL
,0,state
->box
,NULL
);
608 clear_mat(state
->box_rel
);
610 if (PRESERVE_SHAPE(*ir
))
611 do_box_rel(ir
,state
->box_rel
,state
->box
,TRUE
);
614 void preserve_box_shape(t_inputrec
*ir
,matrix box_rel
,matrix b
)
616 if (PRESERVE_SHAPE(*ir
))
617 do_box_rel(ir
,box_rel
,b
,FALSE
);
620 void add_t_atoms(t_atoms
*atoms
,int natom_extra
,int nres_extra
)
626 srenew(atoms
->atomname
,atoms
->nr
+natom_extra
);
627 srenew(atoms
->atom
,atoms
->nr
+natom_extra
);
628 if (NULL
!= atoms
->pdbinfo
)
629 srenew(atoms
->pdbinfo
,atoms
->nr
+natom_extra
);
630 if (NULL
!= atoms
->atomtype
)
631 srenew(atoms
->atomtype
,atoms
->nr
+natom_extra
);
632 if (NULL
!= atoms
->atomtypeB
)
633 srenew(atoms
->atomtypeB
,atoms
->nr
+natom_extra
);
634 for(i
=atoms
->nr
; (i
<atoms
->nr
+natom_extra
); i
++) {
635 atoms
->atomname
[i
] = NULL
;
636 memset(&atoms
->atom
[i
],0,sizeof(atoms
->atom
[i
]));
637 if (NULL
!= atoms
->pdbinfo
)
638 memset(&atoms
->pdbinfo
[i
],0,sizeof(atoms
->pdbinfo
[i
]));
639 if (NULL
!= atoms
->atomtype
)
640 atoms
->atomtype
[i
] = NULL
;
641 if (NULL
!= atoms
->atomtypeB
)
642 atoms
->atomtypeB
[i
] = NULL
;
644 atoms
->nr
+= natom_extra
;
648 srenew(atoms
->resinfo
,atoms
->nres
+nres_extra
);
649 for(i
=atoms
->nres
; (i
<atoms
->nres
+nres_extra
); i
++) {
650 memset(&atoms
->resinfo
[i
],0,sizeof(atoms
->resinfo
[i
]));
652 atoms
->nres
+= nres_extra
;
656 void init_t_atoms(t_atoms
*atoms
, int natoms
, gmx_bool bPdbinfo
)
660 snew(atoms
->atomname
,natoms
);
661 atoms
->atomtype
=NULL
;
662 atoms
->atomtypeB
=NULL
;
663 snew(atoms
->resinfo
,natoms
);
664 snew(atoms
->atom
,natoms
);
666 snew(atoms
->pdbinfo
,natoms
);
671 t_atoms
*copy_t_atoms(t_atoms
*src
)
677 init_t_atoms(dst
,src
->nr
,(NULL
!= src
->pdbinfo
));
679 if (NULL
!= src
->atomname
)
680 snew(dst
->atomname
,src
->nr
);
681 if (NULL
!= src
->atomtype
)
682 snew(dst
->atomtype
,src
->nr
);
683 if (NULL
!= src
->atomtypeB
)
684 snew(dst
->atomtypeB
,src
->nr
);
685 for(i
=0; (i
<src
->nr
); i
++) {
686 dst
->atom
[i
] = src
->atom
[i
];
687 if (NULL
!= src
->pdbinfo
)
688 dst
->pdbinfo
[i
] = src
->pdbinfo
[i
];
689 if (NULL
!= src
->atomname
)
690 dst
->atomname
[i
] = src
->atomname
[i
];
691 if (NULL
!= src
->atomtype
)
692 dst
->atomtype
[i
] = src
->atomtype
[i
];
693 if (NULL
!= src
->atomtypeB
)
694 dst
->atomtypeB
[i
] = src
->atomtypeB
[i
];
696 dst
->nres
= src
->nres
;
697 for(i
=0; (i
<src
->nres
); i
++) {
698 dst
->resinfo
[i
] = src
->resinfo
[i
];
703 void t_atoms_set_resinfo(t_atoms
*atoms
,int atom_ind
,t_symtab
*symtab
,
704 const char *resname
,int resnr
,unsigned char ic
,
705 int chainnum
, char chainid
)
709 ri
= &atoms
->resinfo
[atoms
->atom
[atom_ind
].resind
];
710 ri
->name
= put_symtab(symtab
,resname
);
714 ri
->chainnum
= chainnum
;
715 ri
->chainid
= chainid
;
718 void free_t_atoms(t_atoms
*atoms
,gmx_bool bFreeNames
)
723 for(i
=0; i
<atoms
->nr
; i
++) {
724 sfree(*atoms
->atomname
[i
]);
725 *atoms
->atomname
[i
]=NULL
;
727 for(i
=0; i
<atoms
->nres
; i
++) {
728 sfree(*atoms
->resinfo
[i
].name
);
729 *atoms
->resinfo
[i
].name
=NULL
;
732 sfree(atoms
->atomname
);
733 /* Do we need to free atomtype and atomtypeB as well ? */
734 sfree(atoms
->resinfo
);
737 sfree(atoms
->pdbinfo
);
742 real
max_cutoff(real cutoff1
,real cutoff2
)
744 if (cutoff1
== 0 || cutoff2
== 0)
750 return max(cutoff1
,cutoff2
);