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 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
40 #include <sys/types.h>
64 #include "sortwater.h"
66 #include "gmx_fatal.h"
72 #include "vsite_parm.h"
78 #include "compute_io.h"
79 #include "gpp_atomtype.h"
80 #include "gpp_tomorse.h"
81 #include "mtop_util.h"
84 static int rm_interactions(int ifunc
,int nrmols
,t_molinfo mols
[])
89 /* For all the molecule types */
90 for(i
=0; i
<nrmols
; i
++) {
91 n
+= mols
[i
].plist
[ifunc
].nr
;
92 mols
[i
].plist
[ifunc
].nr
=0;
97 static int check_atom_names(const char *fn1
, const char *fn2
,
98 gmx_mtop_t
*mtop
, t_atoms
*at
)
100 int mb
,m
,i
,j
,nmismatch
;
102 #define MAXMISMATCH 20
104 if (mtop
->natoms
!= at
->nr
)
105 gmx_incons("comparing atom names");
109 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
110 tat
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
111 for(m
=0; m
<mtop
->molblock
[mb
].nmol
; m
++) {
112 for(j
=0; j
< tat
->nr
; j
++) {
113 if (strcmp( *(tat
->atomname
[j
]) , *(at
->atomname
[i
]) ) != 0) {
114 if (nmismatch
< MAXMISMATCH
) {
116 "Warning: atom name %d in %s and %s does not match (%s - %s)\n",
117 i
+1, fn1
, fn2
, *(tat
->atomname
[j
]), *(at
->atomname
[i
]));
118 } else if (nmismatch
== MAXMISMATCH
) {
119 fprintf(stderr
,"(more than %d non-matching atom names)\n",MAXMISMATCH
);
131 static void check_eg_vs_cg(gmx_mtop_t
*mtop
)
133 int astart
,mb
,m
,cg
,j
,firstj
;
134 unsigned char firsteg
,eg
;
137 /* Go through all the charge groups and make sure all their
138 * atoms are in the same energy group.
142 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
143 molt
= &mtop
->moltype
[mtop
->molblock
[mb
].type
];
144 for(m
=0; m
<mtop
->molblock
[mb
].nmol
; m
++) {
145 for(cg
=0; cg
<molt
->cgs
.nr
;cg
++) {
146 /* Get the energy group of the first atom in this charge group */
147 firstj
= astart
+ molt
->cgs
.index
[cg
];
148 firsteg
= ggrpnr(&mtop
->groups
,egcENER
,firstj
);
149 for(j
=molt
->cgs
.index
[cg
]+1;j
<molt
->cgs
.index
[cg
+1];j
++) {
150 eg
= ggrpnr(&mtop
->groups
,egcENER
,astart
+j
);
152 gmx_fatal(FARGS
,"atoms %d and %d in charge group %d of molecule type '%s' are in different energy groups",
153 firstj
+1,astart
+j
+1,cg
+1,*molt
->name
);
157 astart
+= molt
->atoms
.nr
;
162 static void check_cg_sizes(const char *topfn
,t_block
*cgs
,warninp_t wi
)
165 char warn_buf
[STRLEN
];
168 for(cg
=0; cg
<cgs
->nr
; cg
++)
170 maxsize
= max(maxsize
,cgs
->index
[cg
+1]-cgs
->index
[cg
]);
174 gmx_fatal(FARGS
,"The largst charge group contains %d atoms. The maximum is 32.",maxsize
);
175 else if (maxsize
> 10)
177 set_warning_line(wi
,topfn
,-1);
179 "The largest charge group contains %d atoms.\n"
180 "Since atoms only see each other when the centers of geometry of the charge groups they belong to are within the cut-off distance, too large charge groups can lead to serious cut-off artifacts.\n"
181 "For efficiency and accuracy, charge group should consist of a few atoms.\n"
182 "For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc.",
184 warning_note(wi
,warn_buf
);
188 static void check_bonds_timestep(gmx_mtop_t
*mtop
,double dt
,warninp_t wi
)
190 /* This check is not intended to ensure accurate integration,
191 * rather it is to signal mistakes in the mdp settings.
192 * A common mistake is to forget to turn on constraints
193 * for MD after energy minimization with flexible bonds.
194 * This check can also detect too large time steps for flexible water
195 * models, but such errors will often be masked by the constraints
196 * mdp options, which turns flexible water into water with bond constraints,
197 * but without an angle constraint. Unfortunately such incorrect use
198 * of water models can not easily be detected without checking
199 * for specific model names.
201 * The stability limit of leap-frog or velocity verlet is 4.44 steps
202 * per oscillational period.
203 * But accurate bonds distributions are lost far before that limit.
204 * To allow relatively common schemes (although not common with Gromacs)
205 * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
206 * we set the note limit to 10.
208 int min_steps_warn
=5;
209 int min_steps_note
=10;
212 gmx_moltype_t
*moltype
,*w_moltype
;
214 t_ilist
*ilist
,*ilb
,*ilc
,*ils
;
216 int i
,a1
,a2
,w_a1
,w_a2
,j
;
217 real twopi2
,limit2
,fc
,re
,m1
,m2
,period2
,w_period2
;
218 gmx_bool bFound
,bWater
,bWarn
;
219 char warn_buf
[STRLEN
];
221 ip
= mtop
->ffparams
.iparams
;
223 twopi2
= sqr(2*M_PI
);
225 limit2
= sqr(min_steps_note
*dt
);
231 for(molt
=0; molt
<mtop
->nmoltype
; molt
++)
233 moltype
= &mtop
->moltype
[molt
];
234 atom
= moltype
->atoms
.atom
;
235 ilist
= moltype
->ilist
;
236 ilc
= &ilist
[F_CONSTR
];
237 ils
= &ilist
[F_SETTLE
];
238 for(ftype
=0; ftype
<F_NRE
; ftype
++)
240 if (!(ftype
== F_BONDS
|| ftype
== F_G96BONDS
|| ftype
== F_HARMONIC
))
246 for(i
=0; i
<ilb
->nr
; i
+=3)
248 fc
= ip
[ilb
->iatoms
[i
]].harmonic
.krA
;
249 re
= ip
[ilb
->iatoms
[i
]].harmonic
.rA
;
250 if (ftype
== F_G96BONDS
)
252 /* Convert squared sqaure fc to harmonic fc */
255 a1
= ilb
->iatoms
[i
+1];
256 a2
= ilb
->iatoms
[i
+2];
259 if (fc
> 0 && m1
> 0 && m2
> 0)
261 period2
= twopi2
*m1
*m2
/((m1
+ m2
)*fc
);
265 period2
= GMX_FLOAT_MAX
;
269 fprintf(debug
,"fc %g m1 %g m2 %g period %g\n",
270 fc
,m1
,m2
,sqrt(period2
));
272 if (period2
< limit2
)
275 for(j
=0; j
<ilc
->nr
; j
+=3)
277 if ((ilc
->iatoms
[j
+1] == a1
&& ilc
->iatoms
[j
+2] == a2
) ||
278 (ilc
->iatoms
[j
+1] == a2
&& ilc
->iatoms
[j
+2] == a1
))
283 for(j
=0; j
<ils
->nr
; j
+=2)
285 if ((a1
>= ils
->iatoms
[j
+1] && a1
< ils
->iatoms
[j
+1]+3) &&
286 (a2
>= ils
->iatoms
[j
+1] && a2
< ils
->iatoms
[j
+1]+3))
292 (w_moltype
== NULL
|| period2
< w_period2
))
304 if (w_moltype
!= NULL
)
306 bWarn
= (w_period2
< sqr(min_steps_warn
*dt
));
307 /* A check that would recognize most water models */
308 bWater
= ((*w_moltype
->atoms
.atomname
[0])[0] == 'O' &&
309 w_moltype
->atoms
.nr
<= 5);
310 sprintf(warn_buf
,"The bond in molecule-type %s between atoms %d %s and %d %s has an estimated oscillational period of %.1e ps, which is less than %d times the time step of %.1e ps.\n"
313 w_a1
+1,*w_moltype
->atoms
.atomname
[w_a1
],
314 w_a2
+1,*w_moltype
->atoms
.atomname
[w_a2
],
315 sqrt(w_period2
),bWarn
? min_steps_warn
: min_steps_note
,dt
,
317 "Maybe you asked for fexible water." :
318 "Maybe you forgot to change the constraints mdp option.");
321 warning(wi
,warn_buf
);
325 warning_note(wi
,warn_buf
);
330 static void check_vel(gmx_mtop_t
*mtop
,rvec v
[])
332 gmx_mtop_atomloop_all_t aloop
;
336 aloop
= gmx_mtop_atomloop_all_init(mtop
);
337 while (gmx_mtop_atomloop_all_next(aloop
,&a
,&atom
)) {
338 if (atom
->ptype
== eptShell
||
339 atom
->ptype
== eptBond
||
340 atom
->ptype
== eptVSite
) {
346 static gmx_bool
nint_ftype(gmx_mtop_t
*mtop
,t_molinfo
*mi
,int ftype
)
351 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
352 nint
+= mtop
->molblock
[mb
].nmol
*mi
[mtop
->molblock
[mb
].type
].plist
[ftype
].nr
;
358 /* This routine reorders the molecule type array
359 * in the order of use in the molblocks,
360 * unused molecule types are deleted.
362 static void renumber_moltypes(gmx_mtop_t
*sys
,
363 int *nmolinfo
,t_molinfo
**molinfo
)
369 snew(order
,*nmolinfo
);
371 for(mb
=0; mb
<sys
->nmolblock
; mb
++) {
372 for(i
=0; i
<norder
; i
++) {
373 if (order
[i
] == sys
->molblock
[mb
].type
) {
378 /* This type did not occur yet, add it */
379 order
[norder
] = sys
->molblock
[mb
].type
;
380 /* Renumber the moltype in the topology */
383 sys
->molblock
[mb
].type
= i
;
386 /* We still need to reorder the molinfo structs */
388 for(mi
=0; mi
<*nmolinfo
; mi
++) {
389 for(i
=0; i
<norder
; i
++) {
390 if (order
[i
] == mi
) {
395 done_mi(&(*molinfo
)[mi
]);
397 minew
[i
] = (*molinfo
)[mi
];
406 static void molinfo2mtop(int nmi
,t_molinfo
*mi
,gmx_mtop_t
*mtop
)
411 mtop
->nmoltype
= nmi
;
412 snew(mtop
->moltype
,nmi
);
413 for(m
=0; m
<nmi
; m
++) {
414 molt
= &mtop
->moltype
[m
];
415 molt
->name
= mi
[m
].name
;
416 molt
->atoms
= mi
[m
].atoms
;
417 /* ilists are copied later */
418 molt
->cgs
= mi
[m
].cgs
;
419 molt
->excls
= mi
[m
].excls
;
424 new_status(const char *topfile
,const char *topppfile
,const char *confin
,
425 t_gromppopts
*opts
,t_inputrec
*ir
,gmx_bool bZero
,
426 gmx_bool bGenVel
,gmx_bool bVerbose
,t_state
*state
,
427 gpp_atomtype_t atype
,gmx_mtop_t
*sys
,
428 int *nmi
,t_molinfo
**mi
,t_params plist
[],
429 int *comb
,double *reppow
,real
*fudgeQQ
,
433 t_molinfo
*molinfo
=NULL
;
435 gmx_molblock_t
*molblock
,*molbs
;
437 int mb
,i
,nrmols
,nmismatch
;
440 char warn_buf
[STRLEN
];
444 /* Set gmx_boolean for GB */
445 if(ir
->implicit_solvent
)
448 /* TOPOLOGY processing */
449 sys
->name
= do_top(bVerbose
,topfile
,topppfile
,opts
,bZero
,&(sys
->symtab
),
450 plist
,comb
,reppow
,fudgeQQ
,
451 atype
,&nrmols
,&molinfo
,ir
,
452 &nmolblock
,&molblock
,bGB
,
456 snew(sys
->molblock
,nmolblock
);
459 for(mb
=0; mb
<nmolblock
; mb
++) {
460 if (sys
->nmolblock
> 0 &&
461 molblock
[mb
].type
== sys
->molblock
[sys
->nmolblock
-1].type
) {
462 /* Merge consecutive blocks with the same molecule type */
463 sys
->molblock
[sys
->nmolblock
-1].nmol
+= molblock
[mb
].nmol
;
464 sys
->natoms
+= molblock
[mb
].nmol
*sys
->molblock
[sys
->nmolblock
-1].natoms_mol
;
465 } else if (molblock
[mb
].nmol
> 0) {
466 /* Add a new molblock to the topology */
467 molbs
= &sys
->molblock
[sys
->nmolblock
];
468 *molbs
= molblock
[mb
];
469 molbs
->natoms_mol
= molinfo
[molbs
->type
].atoms
.nr
;
470 molbs
->nposres_xA
= 0;
471 molbs
->nposres_xB
= 0;
472 sys
->natoms
+= molbs
->nmol
*molbs
->natoms_mol
;
476 if (sys
->nmolblock
== 0) {
477 gmx_fatal(FARGS
,"No molecules were defined in the system");
480 renumber_moltypes(sys
,&nrmols
,&molinfo
);
483 convert_harmonics(nrmols
,molinfo
,atype
);
485 if (ir
->eDisre
== edrNone
) {
486 i
= rm_interactions(F_DISRES
,nrmols
,molinfo
);
488 set_warning_line(wi
,"unknown",-1);
489 sprintf(warn_buf
,"disre = no, removed %d distance restraints",i
);
490 warning_note(wi
,warn_buf
);
493 if (opts
->bOrire
== FALSE
) {
494 i
= rm_interactions(F_ORIRES
,nrmols
,molinfo
);
496 set_warning_line(wi
,"unknown",-1);
497 sprintf(warn_buf
,"orire = no, removed %d orientation restraints",i
);
498 warning_note(wi
,warn_buf
);
501 if (opts
->bDihre
== FALSE
) {
502 i
= rm_interactions(F_DIHRES
,nrmols
,molinfo
);
504 set_warning_line(wi
,"unknown",-1);
505 sprintf(warn_buf
,"dihre = no, removed %d dihedral restraints",i
);
506 warning_note(wi
,warn_buf
);
510 /* Copy structures from msys to sys */
511 molinfo2mtop(nrmols
,molinfo
,sys
);
513 gmx_mtop_finalize(sys
);
515 /* COORDINATE file processing */
517 fprintf(stderr
,"processing coordinates...\n");
519 get_stx_coordnum(confin
,&state
->natoms
);
520 if (state
->natoms
!= sys
->natoms
)
521 gmx_fatal(FARGS
,"number of coordinates in coordinate file (%s, %d)\n"
522 " does not match topology (%s, %d)",
523 confin
,state
->natoms
,topfile
,sys
->natoms
);
525 /* make space for coordinates and velocities */
528 init_t_atoms(confat
,state
->natoms
,FALSE
);
529 init_state(state
,state
->natoms
,0,0,0);
530 read_stx_conf(confin
,title
,confat
,state
->x
,state
->v
,NULL
,state
->box
);
531 /* This call fixes the box shape for runs with pressure scaling */
532 set_box_rel(ir
,state
);
534 nmismatch
= check_atom_names(topfile
, confin
, sys
, confat
);
535 free_t_atoms(confat
,TRUE
);
539 sprintf(buf
,"%d non-matching atom name%s\n"
540 "atom names from %s will be used\n"
541 "atom names from %s will be ignored\n",
542 nmismatch
,(nmismatch
== 1) ? "" : "s",topfile
,confin
);
546 fprintf(stderr
,"double-checking input for internal consistency...\n");
547 double_check(ir
,state
->box
,nint_ftype(sys
,molinfo
,F_CONSTR
),wi
);
552 gmx_mtop_atomloop_all_t aloop
;
555 snew(mass
,state
->natoms
);
556 aloop
= gmx_mtop_atomloop_all_init(sys
);
557 while (gmx_mtop_atomloop_all_next(aloop
,&i
,&atom
)) {
561 if (opts
->seed
== -1) {
562 opts
->seed
= make_seed();
563 fprintf(stderr
,"Setting gen_seed to %d\n",opts
->seed
);
565 maxwell_speed(opts
->tempi
,opts
->seed
,sys
,state
->v
);
567 stop_cm(stdout
,state
->natoms
,mass
,state
->x
,state
->v
);
575 static void copy_state(const char *slog
,t_trxframe
*fr
,
576 gmx_bool bReadVel
,t_state
*state
,
581 if (fr
->not_ok
& FRAME_NOT_OK
)
583 gmx_fatal(FARGS
,"Can not start from an incomplete frame");
587 gmx_fatal(FARGS
,"Did not find a frame with coordinates in file %s",
591 for(i
=0; i
<state
->natoms
; i
++)
593 copy_rvec(fr
->x
[i
],state
->x
[i
]);
599 gmx_incons("Trajecory frame unexpectedly does not contain velocities");
601 for(i
=0; i
<state
->natoms
; i
++)
603 copy_rvec(fr
->v
[i
],state
->v
[i
]);
608 copy_mat(fr
->box
,state
->box
);
611 *use_time
= fr
->time
;
614 static void cont_status(const char *slog
,const char *ener
,
615 gmx_bool bNeedVel
,gmx_bool bGenVel
, real fr_time
,
616 t_inputrec
*ir
,t_state
*state
,
618 const output_env_t oenv
)
619 /* If fr_time == -1 read the last frame available which is complete */
627 bReadVel
= (bNeedVel
&& !bGenVel
);
630 "Reading Coordinates%s and Box size from old trajectory\n",
631 bReadVel
? ", Velocities" : "");
634 fprintf(stderr
,"Will read whole trajectory\n");
638 fprintf(stderr
,"Will read till time %g\n",fr_time
);
644 fprintf(stderr
,"Velocities generated: "
645 "ignoring velocities in input trajectory\n");
647 read_first_frame(oenv
,&fp
,slog
,&fr
,TRX_NEED_X
);
651 read_first_frame(oenv
,&fp
,slog
,&fr
,TRX_NEED_X
| TRX_NEED_V
);
657 "WARNING: Did not find a frame with velocities in file %s,\n"
658 " all velocities will be set to zero!\n\n",slog
);
659 for(i
=0; i
<sys
->natoms
; i
++)
661 clear_rvec(state
->v
[i
]);
664 /* Search for a frame without velocities */
666 read_first_frame(oenv
,&fp
,slog
,&fr
,TRX_NEED_X
);
670 state
->natoms
= fr
.natoms
;
672 if (sys
->natoms
!= state
->natoms
)
674 gmx_fatal(FARGS
,"Number of atoms in Topology "
675 "is not the same as in Trajectory");
677 copy_state(slog
,&fr
,bReadVel
,state
,&use_time
);
679 /* Find the appropriate frame */
680 while ((fr_time
== -1 || fr
.time
< fr_time
) &&
681 read_next_frame(oenv
,fp
,&fr
))
683 copy_state(slog
,&fr
,bReadVel
,state
,&use_time
);
688 /* Set the relative box lengths for preserving the box shape.
689 * Note that this call can lead to differences in the last bit
690 * with respect to using tpbconv to create a tpx file.
692 set_box_rel(ir
,state
);
694 fprintf(stderr
,"Using frame at t = %g ps\n",use_time
);
695 fprintf(stderr
,"Starting time for run is %g ps\n",ir
->init_t
);
697 if ((ir
->epc
!= epcNO
|| ir
->etc
==etcNOSEHOOVER
) && ener
)
699 get_enx_state(ener
,use_time
,&sys
->groups
,ir
,state
);
700 preserve_box_shape(ir
,state
->box_rel
,state
->boxv
);
704 static void read_posres(gmx_mtop_t
*mtop
,t_molinfo
*molinfo
,gmx_bool bTopB
,
706 int rc_scaling
, int ePBC
,
710 gmx_bool bFirst
= TRUE
;
716 int natoms
,npbcdim
=0;
717 char warn_buf
[STRLEN
],title
[STRLEN
];
718 int a
,i
,ai
,j
,k
,mb
,nat_molb
;
719 gmx_molblock_t
*molb
;
723 get_stx_coordnum(fn
,&natoms
);
724 if (natoms
!= mtop
->natoms
) {
725 sprintf(warn_buf
,"The number of atoms in %s (%d) does not match the number of atoms in the topology (%d). Will assume that the first %d atoms in the topology and %s match.",fn
,natoms
,mtop
->natoms
,min(mtop
->natoms
,natoms
),fn
);
726 warning(wi
,warn_buf
);
730 init_t_atoms(&dumat
,natoms
,FALSE
);
731 read_stx_conf(fn
,title
,&dumat
,x
,v
,NULL
,box
);
733 npbcdim
= ePBC2npbcdim(ePBC
);
735 if (rc_scaling
!= erscNO
) {
736 copy_mat(box
,invbox
);
737 for(j
=npbcdim
; j
<DIM
; j
++) {
738 clear_rvec(invbox
[j
]);
741 m_inv_ur0(invbox
,invbox
);
744 /* Copy the reference coordinates to mtop */
748 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
749 molb
= &mtop
->molblock
[mb
];
750 nat_molb
= molb
->nmol
*mtop
->moltype
[molb
->type
].atoms
.nr
;
751 pr
= &(molinfo
[molb
->type
].plist
[F_POSRES
]);
753 atom
= mtop
->moltype
[molb
->type
].atoms
.atom
;
754 for(i
=0; (i
<pr
->nr
); i
++) {
757 gmx_fatal(FARGS
,"Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
758 ai
+1,*molinfo
[molb
->type
].name
,fn
,natoms
);
760 if (rc_scaling
== erscCOM
) {
761 /* Determine the center of mass of the posres reference coordinates */
762 for(j
=0; j
<npbcdim
; j
++) {
763 sum
[j
] += atom
[ai
].m
*x
[a
+ai
][j
];
765 totmass
+= atom
[ai
].m
;
769 molb
->nposres_xA
= nat_molb
;
770 snew(molb
->posres_xA
,molb
->nposres_xA
);
771 for(i
=0; i
<nat_molb
; i
++) {
772 copy_rvec(x
[a
+i
],molb
->posres_xA
[i
]);
775 molb
->nposres_xB
= nat_molb
;
776 snew(molb
->posres_xB
,molb
->nposres_xB
);
777 for(i
=0; i
<nat_molb
; i
++) {
778 copy_rvec(x
[a
+i
],molb
->posres_xB
[i
]);
784 if (rc_scaling
== erscCOM
) {
786 gmx_fatal(FARGS
,"The total mass of the position restraint atoms is 0");
787 for(j
=0; j
<npbcdim
; j
++)
788 com
[j
] = sum
[j
]/totmass
;
789 fprintf(stderr
,"The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f\n",com
[XX
],com
[YY
],com
[ZZ
]);
792 if (rc_scaling
!= erscNO
) {
793 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
794 molb
= &mtop
->molblock
[mb
];
795 nat_molb
= molb
->nmol
*mtop
->moltype
[molb
->type
].atoms
.nr
;
796 if (molb
->nposres_xA
> 0 || molb
->nposres_xB
> 0) {
797 xp
= (!bTopB
? molb
->posres_xA
: molb
->posres_xB
);
798 for(i
=0; i
<nat_molb
; i
++) {
799 for(j
=0; j
<npbcdim
; j
++) {
800 if (rc_scaling
== erscALL
) {
801 /* Convert from Cartesian to crystal coordinates */
802 xp
[i
][j
] *= invbox
[j
][j
];
803 for(k
=j
+1; k
<npbcdim
; k
++) {
804 xp
[i
][j
] += invbox
[k
][j
]*xp
[i
][k
];
806 } else if (rc_scaling
== erscCOM
) {
807 /* Subtract the center of mass */
815 if (rc_scaling
== erscCOM
) {
816 /* Convert the COM from Cartesian to crystal coordinates */
817 for(j
=0; j
<npbcdim
; j
++) {
818 com
[j
] *= invbox
[j
][j
];
819 for(k
=j
+1; k
<npbcdim
; k
++) {
820 com
[j
] += invbox
[k
][j
]*com
[k
];
826 free_t_atoms(&dumat
,TRUE
);
831 static void gen_posres(gmx_mtop_t
*mtop
,t_molinfo
*mi
,
832 char *fnA
, char *fnB
,
833 int rc_scaling
, int ePBC
,
839 read_posres (mtop
,mi
,FALSE
,fnA
,rc_scaling
,ePBC
,com
,wi
);
840 if (strcmp(fnA
,fnB
) != 0) {
841 read_posres(mtop
,mi
,TRUE
,fnB
,rc_scaling
,ePBC
,comB
,wi
);
845 static void set_wall_atomtype(gpp_atomtype_t at
,t_gromppopts
*opts
,
851 fprintf(stderr
,"Searching the wall atom type(s)\n");
852 for(i
=0; i
<ir
->nwall
; i
++)
853 ir
->wall_atomtype
[i
] = get_atomtype_type(opts
->wall_atomtype
[i
],at
);
856 static int nrdf_internal(t_atoms
*atoms
)
861 for(i
=0; i
<atoms
->nr
; i
++) {
862 /* Vsite ptype might not be set here yet, so also check the mass */
863 if ((atoms
->atom
[i
].ptype
== eptAtom
||
864 atoms
->atom
[i
].ptype
== eptNucleus
)
865 && atoms
->atom
[i
].m
> 0) {
870 case 0: nrdf
= 0; break;
871 case 1: nrdf
= 0; break;
872 case 2: nrdf
= 1; break;
873 default: nrdf
= nmass
*3 - 6; break;
896 q
= (y
[i
+1]-2.0*y
[i
]+y
[i
-1])/dx
;
897 u
[i
] = (3.0*q
/dx
-0.5*u
[i
-1])/p
;
904 y2
[i
] = y2
[i
]*y2
[i
+1]+u
[i
];
910 interpolate1d( double xmin
,
923 a
= (xmin
+(ix
+1)*dx
-x
)/dx
;
924 b
= (x
-xmin
-ix
*dx
)/dx
;
926 *y
= a
*ya
[ix
]+b
*ya
[ix
+1]+((a
*a
*a
-a
)*y2a
[ix
]+(b
*b
*b
-b
)*y2a
[ix
+1])*(dx
*dx
)/6.0;
927 *y1
= (ya
[ix
+1]-ya
[ix
])/dx
-(3.0*a
*a
-1.0)/6.0*dx
*y2a
[ix
]+(3.0*b
*b
-1.0)/6.0*dx
*y2a
[ix
+1];
932 setup_cmap (int grid_spacing
,
935 gmx_cmap_t
* cmap_grid
)
937 double *tmp_u
,*tmp_u2
,*tmp_yy
,*tmp_y1
,*tmp_t2
,*tmp_grid
;
939 int i
,j
,k
,ii
,jj
,kk
,idx
;
941 double dx
,xmin
,v
,v1
,v2
,v12
;
944 snew(tmp_u
,2*grid_spacing
);
945 snew(tmp_u2
,2*grid_spacing
);
946 snew(tmp_yy
,2*grid_spacing
);
947 snew(tmp_y1
,2*grid_spacing
);
948 snew(tmp_t2
,2*grid_spacing
*2*grid_spacing
);
949 snew(tmp_grid
,2*grid_spacing
*2*grid_spacing
);
951 dx
= 360.0/grid_spacing
;
952 xmin
= -180.0-dx
*grid_spacing
/2;
956 /* Compute an offset depending on which cmap we are using
957 * Offset will be the map number multiplied with the grid_spacing * grid_spacing * 2
959 offset
= kk
* grid_spacing
* grid_spacing
* 2;
961 for(i
=0;i
<2*grid_spacing
;i
++)
963 ii
=(i
+grid_spacing
-grid_spacing
/2)%grid_spacing
;
965 for(j
=0;j
<2*grid_spacing
;j
++)
967 jj
=(j
+grid_spacing
-grid_spacing
/2)%grid_spacing
;
968 tmp_grid
[i
*grid_spacing
*2+j
] = grid
[offset
+ii
*grid_spacing
+jj
];
972 for(i
=0;i
<2*grid_spacing
;i
++)
974 spline1d(dx
,&(tmp_grid
[2*grid_spacing
*i
]),2*grid_spacing
,tmp_u
,&(tmp_t2
[2*grid_spacing
*i
]));
977 for(i
=grid_spacing
/2;i
<grid_spacing
+grid_spacing
/2;i
++)
979 ii
= i
-grid_spacing
/2;
982 for(j
=grid_spacing
/2;j
<grid_spacing
+grid_spacing
/2;j
++)
984 jj
= j
-grid_spacing
/2;
987 for(k
=0;k
<2*grid_spacing
;k
++)
989 interpolate1d(xmin
,dx
,&(tmp_grid
[2*grid_spacing
*k
]),
990 &(tmp_t2
[2*grid_spacing
*k
]),psi
,&tmp_yy
[k
],&tmp_y1
[k
]);
993 spline1d(dx
,tmp_yy
,2*grid_spacing
,tmp_u
,tmp_u2
);
994 interpolate1d(xmin
,dx
,tmp_yy
,tmp_u2
,phi
,&v
,&v1
);
995 spline1d(dx
,tmp_y1
,2*grid_spacing
,tmp_u
,tmp_u2
);
996 interpolate1d(xmin
,dx
,tmp_y1
,tmp_u2
,phi
,&v2
,&v12
);
998 idx
= ii
*grid_spacing
+jj
;
999 cmap_grid
->cmapdata
[kk
].cmap
[idx
*4] = grid
[offset
+ii
*grid_spacing
+jj
];
1000 cmap_grid
->cmapdata
[kk
].cmap
[idx
*4+1] = v1
;
1001 cmap_grid
->cmapdata
[kk
].cmap
[idx
*4+2] = v2
;
1002 cmap_grid
->cmapdata
[kk
].cmap
[idx
*4+3] = v12
;
1008 void init_cmap_grid(gmx_cmap_t
*cmap_grid
, int ngrid
, int grid_spacing
)
1012 cmap_grid
->ngrid
= ngrid
;
1013 cmap_grid
->grid_spacing
= grid_spacing
;
1014 nelem
= cmap_grid
->grid_spacing
*cmap_grid
->grid_spacing
;
1016 snew(cmap_grid
->cmapdata
,ngrid
);
1018 for(i
=0;i
<cmap_grid
->ngrid
;i
++)
1020 snew(cmap_grid
->cmapdata
[i
].cmap
,4*nelem
);
1025 static int count_constraints(gmx_mtop_t
*mtop
,t_molinfo
*mi
,warninp_t wi
)
1027 int count
,count_mol
,i
,mb
;
1028 gmx_molblock_t
*molb
;
1033 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
1035 molb
= &mtop
->molblock
[mb
];
1036 plist
= mi
[molb
->type
].plist
;
1038 for(i
=0; i
<F_NRE
; i
++) {
1040 count_mol
+= 3*plist
[i
].nr
;
1041 else if (interaction_function
[i
].flags
& IF_CONSTRAINT
)
1042 count_mol
+= plist
[i
].nr
;
1045 if (count_mol
> nrdf_internal(&mi
[molb
->type
].atoms
)) {
1047 "Molecule type '%s' has %d constraints.\n"
1048 "For stability and efficiency there should not be more constraints than internal number of degrees of freedom: %d.\n",
1049 *mi
[molb
->type
].name
,count_mol
,
1050 nrdf_internal(&mi
[molb
->type
].atoms
));
1053 count
+= molb
->nmol
*count_mol
;
1059 static void check_gbsa_charges_lj(gmx_mtop_t
*sys
, gpp_atomtype_t atype
)
1061 int i
,nmiss
,natoms
,mt
;
1063 const t_atoms
*atoms
;
1067 for(mt
=0;mt
<sys
->nmoltype
;mt
++)
1069 atoms
= &sys
->moltype
[mt
].atoms
;
1072 for(i
=0;i
<natoms
;i
++)
1075 q
= atoms
->atom
[i
].q
;
1076 if( (get_atomtype_radius(atoms
->atom
[i
].type
,atype
) <= 0 ||
1077 get_atomtype_vol(atoms
->atom
[i
].type
,atype
) <= 0 ||
1078 get_atomtype_surftens(atoms
->atom
[i
].type
,atype
) <= 0 ||
1079 get_atomtype_gb_radius(atoms
->atom
[i
].type
,atype
) <= 0 ||
1080 get_atomtype_S_hct(atoms
->atom
[i
].type
,atype
) <= 0) && q
!=0)
1082 fprintf(stderr
,"GB parameter(s) missing or negative for atom type '%s' while charge is %g\n",
1083 get_atomtype_name(atoms
->atom
[i
].type
,atype
),q
);
1091 gmx_fatal(FARGS
,"Can't do GB electrostatics; the forcefield is missing %d values for\n"
1092 "particles that also have a charge\n.",nmiss
);
1098 static void check_gbsa_params(t_inputrec
*ir
,gpp_atomtype_t atype
)
1102 /* If we are doing GBSA, check that we got the parameters we need
1103 * This checking is to see if there are GBSA paratmeters for all
1104 * atoms in the force field. To go around this for testing purposes
1105 * comment out the nerror++ counter temporarily
1108 for(i
=0;i
<get_atomtype_ntypes(atype
);i
++)
1110 if (get_atomtype_radius(i
,atype
) < 0 ||
1111 get_atomtype_vol(i
,atype
) < 0 ||
1112 get_atomtype_surftens(i
,atype
) < 0 ||
1113 get_atomtype_gb_radius(i
,atype
) < 0 ||
1114 get_atomtype_S_hct(i
,atype
) < 0)
1116 fprintf(stderr
,"GB parameter(s) missing or negative for atom type '%s'\n",
1117 get_atomtype_name(i
,atype
));
1124 gmx_fatal(FARGS
,"Can't do GB electrostatics; the forcefield is missing %d values for\n"
1125 "atomtype radii, or they might be negative\n.",nmiss
);
1130 int main (int argc
, char *argv
[])
1132 static const char *desc
[] = {
1133 "The gromacs preprocessor",
1134 "reads a molecular topology file, checks the validity of the",
1135 "file, expands the topology from a molecular description to an atomic",
1136 "description. The topology file contains information about",
1137 "molecule types and the number of molecules, the preprocessor",
1138 "copies each molecule as needed. ",
1139 "There is no limitation on the number of molecule types. ",
1140 "Bonds and bond-angles can be converted into constraints, separately",
1141 "for hydrogens and heavy atoms.",
1142 "Then a coordinate file is read and velocities can be generated",
1143 "from a Maxwellian distribution if requested.",
1144 "grompp also reads parameters for the mdrun ",
1145 "(eg. number of MD steps, time step, cut-off), and others such as",
1146 "NEMD parameters, which are corrected so that the net acceleration",
1148 "Eventually a binary file is produced that can serve as the sole input",
1149 "file for the MD program.[PAR]",
1151 "grompp uses the atom names from the topology file. The atom names",
1152 "in the coordinate file (option [TT]-c[tt]) are only read to generate",
1153 "warnings when they do not match the atom names in the topology.",
1154 "Note that the atom names are irrelevant for the simulation as",
1155 "only the atom types are used for generating interaction parameters.[PAR]",
1157 "grompp uses a built-in preprocessor to resolve includes, macros ",
1158 "etcetera. The preprocessor supports the following keywords:[BR]",
1159 "#ifdef VARIABLE[BR]",
1160 "#ifndef VARIABLE[BR]",
1163 "#define VARIABLE[BR]",
1164 "#undef VARIABLE[BR]"
1165 "#include \"filename\"[BR]",
1166 "#include <filename>[BR]",
1167 "The functioning of these statements in your topology may be modulated by",
1168 "using the following two flags in your [TT]mdp[tt] file:[BR]",
1169 "define = -DVARIABLE1 -DVARIABLE2[BR]",
1170 "include = -I/home/john/doe[BR]",
1171 "For further information a C-programming textbook may help you out.",
1172 "Specifying the [TT]-pp[tt] flag will get the pre-processed",
1173 "topology file written out so that you can verify its contents.[PAR]",
1175 "If your system does not have a c-preprocessor, you can still",
1176 "use grompp, but you do not have access to the features ",
1177 "from the cpp. Command line options to the c-preprocessor can be given",
1178 "in the [TT].mdp[tt] file. See your local manual (man cpp).[PAR]",
1180 "When using position restraints a file with restraint coordinates",
1181 "can be supplied with [TT]-r[tt], otherwise restraining will be done",
1182 "with respect to the conformation from the [TT]-c[tt] option.",
1183 "For free energy calculation the the coordinates for the B topology",
1184 "can be supplied with [TT]-rb[tt], otherwise they will be equal to",
1185 "those of the A topology.[PAR]",
1187 "Starting coordinates can be read from trajectory with [TT]-t[tt].",
1188 "The last frame with coordinates and velocities will be read,",
1189 "unless the [TT]-time[tt] option is used.",
1190 "Note that these velocities will not be used when [TT]gen_vel = yes[tt]",
1191 "in your [TT].mdp[tt] file. An energy file can be supplied with",
1192 "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling",
1193 "variables. Note that for continuation it is better and easier to supply",
1194 "a checkpoint file directly to mdrun, since that always contains",
1195 "the complete state of the system and you don't need to generate",
1196 "a new run input file. Note that if you only want to change the number",
1197 "of run steps tpbconv is more convenient than grompp.[PAR]",
1199 "By default all bonded interactions which have constant energy due to",
1200 "virtual site constructions will be removed. If this constant energy is",
1201 "not zero, this will result in a shift in the total energy. All bonded",
1202 "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,",
1203 "all constraints for distances which will be constant anyway because",
1204 "of virtual site constructions will be removed. If any constraints remain",
1205 "which involve virtual sites, a fatal error will result.[PAR]"
1207 "To verify your run input file, please make notice of all warnings",
1208 "on the screen, and correct where necessary. Do also look at the contents",
1209 "of the [TT]mdout.mdp[tt] file, this contains comment lines, as well as",
1210 "the input that [TT]grompp[tt] has read. If in doubt you can start grompp",
1211 "with the [TT]-debug[tt] option which will give you more information",
1212 "in a file called grompp.log (along with real debug info). Finally, you",
1213 "can see the contents of the run input file with the [TT]gmxdump[tt]",
1220 gpp_atomtype_t atype
;
1222 int natoms
,nvsite
,comb
,mt
;
1226 real max_spacing
,fudgeQQ
;
1228 char fn
[STRLEN
],fnB
[STRLEN
];
1229 const char *mdparin
;
1231 gmx_bool bNeedVel
,bGenVel
;
1232 gmx_bool have_atomnumber
;
1234 t_params
*gb_plist
= NULL
;
1235 gmx_genborn_t
*born
= NULL
;
1237 gmx_bool bVerbose
= FALSE
;
1239 char warn_buf
[STRLEN
];
1242 { efMDP
, NULL
, NULL
, ffREAD
},
1243 { efMDP
, "-po", "mdout", ffWRITE
},
1244 { efSTX
, "-c", NULL
, ffREAD
},
1245 { efSTX
, "-r", NULL
, ffOPTRD
},
1246 { efSTX
, "-rb", NULL
, ffOPTRD
},
1247 { efNDX
, NULL
, NULL
, ffOPTRD
},
1248 { efTOP
, NULL
, NULL
, ffREAD
},
1249 { efTOP
, "-pp", "processed", ffOPTWR
},
1250 { efTPX
, "-o", NULL
, ffWRITE
},
1251 { efTRN
, "-t", NULL
, ffOPTRD
},
1252 { efEDR
, "-e", NULL
, ffOPTRD
}
1254 #define NFILE asize(fnm)
1256 /* Command line options */
1257 static gmx_bool bRenum
=TRUE
;
1258 static gmx_bool bRmVSBds
=TRUE
,bZero
=FALSE
;
1259 static int i
,maxwarn
=0;
1260 static real fr_time
=-1;
1262 { "-v", FALSE
, etBOOL
,{&bVerbose
},
1263 "Be loud and noisy" },
1264 { "-time", FALSE
, etREAL
, {&fr_time
},
1265 "Take frame at or first after this time." },
1266 { "-rmvsbds",FALSE
, etBOOL
, {&bRmVSBds
},
1267 "Remove constant bonded interactions with virtual sites" },
1268 { "-maxwarn", FALSE
, etINT
, {&maxwarn
},
1269 "Number of allowed warnings during input processing" },
1270 { "-zero", FALSE
, etBOOL
, {&bZero
},
1271 "Set parameters for bonded interactions without defaults to zero instead of generating an error" },
1272 { "-renum", FALSE
, etBOOL
, {&bRenum
},
1273 "Renumber atomtypes and minimize number of atomtypes" }
1276 CopyRight(stdout
,argv
[0]);
1278 /* Initiate some variables */
1283 /* Parse the command line */
1284 parse_common_args(&argc
,argv
,0,NFILE
,fnm
,asize(pa
),pa
,
1285 asize(desc
),desc
,0,NULL
,&oenv
);
1287 wi
= init_warning(TRUE
,maxwarn
);
1289 /* PARAMETER file processing */
1290 mdparin
= opt2fn("-f",NFILE
,fnm
);
1291 set_warning_line(wi
,mdparin
,-1);
1292 get_ir(mdparin
,opt2fn("-po",NFILE
,fnm
),ir
,opts
,wi
);
1295 fprintf(stderr
,"checking input for internal consistency...\n");
1296 check_ir(mdparin
,ir
,opts
,wi
);
1298 if (ir
->ld_seed
== -1) {
1299 ir
->ld_seed
= make_seed();
1300 fprintf(stderr
,"Setting the LD random seed to %d\n",ir
->ld_seed
);
1303 bNeedVel
= EI_STATE_VELOCITY(ir
->eI
);
1304 bGenVel
= (bNeedVel
&& opts
->bGenVel
);
1309 atype
= init_atomtype();
1311 pr_symtab(debug
,0,"Just opened",&sys
->symtab
);
1313 strcpy(fn
,ftp2fn(efTOP
,NFILE
,fnm
));
1314 if (!gmx_fexist(fn
))
1315 gmx_fatal(FARGS
,"%s does not exist",fn
);
1316 new_status(fn
,opt2fn_null("-pp",NFILE
,fnm
),opt2fn("-c",NFILE
,fnm
),
1317 opts
,ir
,bZero
,bGenVel
,bVerbose
,&state
,
1318 atype
,sys
,&nmi
,&mi
,plist
,&comb
,&reppow
,&fudgeQQ
,
1323 pr_symtab(debug
,0,"After new_status",&sys
->symtab
);
1325 if (count_constraints(sys
,mi
,wi
) && (ir
->eConstrAlg
== econtSHAKE
)) {
1326 if (ir
->eI
== eiCG
|| ir
->eI
== eiLBFGS
) {
1327 sprintf(warn_buf
,"Can not do %s with %s, use %s",
1328 EI(ir
->eI
),econstr_names
[econtSHAKE
],econstr_names
[econtLINCS
]);
1329 warning_error(wi
,warn_buf
);
1331 if (ir
->bPeriodicMols
) {
1332 sprintf(warn_buf
,"Can not do periodic molecules with %s, use %s",
1333 econstr_names
[econtSHAKE
],econstr_names
[econtLINCS
]);
1334 warning_error(wi
,warn_buf
);
1338 /* If we are doing QM/MM, check that we got the atom numbers */
1339 have_atomnumber
= TRUE
;
1340 for (i
=0; i
<get_atomtype_ntypes(atype
); i
++) {
1341 have_atomnumber
= have_atomnumber
&& (get_atomtype_atomnumber(i
,atype
) >= 0);
1343 if (!have_atomnumber
&& ir
->bQMMM
)
1347 "It appears as if you are trying to run a QM/MM calculation, but the force\n"
1348 "field you are using does not contain atom numbers fields. This is an\n"
1349 "optional field (introduced in Gromacs 3.3) for general runs, but mandatory\n"
1350 "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
1351 "an integer just before the mass column in ffXXXnb.itp.\n"
1352 "NB: United atoms have the same atom numbers as normal ones.\n\n");
1355 /* Check for errors in the input now, since they might cause problems
1356 * during processing further down.
1358 check_warning_error(wi
,FARGS
);
1360 if (opt2bSet("-r",NFILE
,fnm
))
1361 sprintf(fn
,"%s",opt2fn("-r",NFILE
,fnm
));
1363 sprintf(fn
,"%s",opt2fn("-c",NFILE
,fnm
));
1364 if (opt2bSet("-rb",NFILE
,fnm
))
1365 sprintf(fnB
,"%s",opt2fn("-rb",NFILE
,fnm
));
1369 if (nint_ftype(sys
,mi
,F_POSRES
) > 0)
1373 fprintf(stderr
,"Reading position restraint coords from %s",fn
);
1374 if (strcmp(fn
,fnB
) == 0)
1376 fprintf(stderr
,"\n");
1380 fprintf(stderr
," and %s\n",fnB
);
1381 if (ir
->efep
!= efepNO
&& ir
->n_flambda
> 0)
1383 warning_error(wi
,"Can not change the position restraint reference coordinates with lambda togther with foreign lambda calculation.");
1387 gen_posres(sys
,mi
,fn
,fnB
,
1388 ir
->refcoord_scaling
,ir
->ePBC
,
1389 ir
->posres_com
,ir
->posres_comB
,
1394 /* set parameters for virtual site construction (not for vsiten) */
1395 for(mt
=0; mt
<sys
->nmoltype
; mt
++) {
1397 set_vsites(bVerbose
, &sys
->moltype
[mt
].atoms
, atype
, mi
[mt
].plist
);
1399 /* now throw away all obsolete bonds, angles and dihedrals: */
1400 /* note: constraints are ALWAYS removed */
1402 for(mt
=0; mt
<sys
->nmoltype
; mt
++) {
1403 clean_vsite_bondeds(mi
[mt
].plist
,sys
->moltype
[mt
].atoms
.nr
,bRmVSBds
);
1407 /* If we are using CMAP, setup the pre-interpolation grid */
1410 init_cmap_grid(&sys
->ffparams
.cmap_grid
, plist
->nc
, plist
->grid_spacing
);
1411 setup_cmap(plist
->grid_spacing
, plist
->nc
, plist
->cmap
,&sys
->ffparams
.cmap_grid
);
1414 set_wall_atomtype(atype
,opts
,ir
);
1416 renum_atype(plist
, sys
, ir
->wall_atomtype
, atype
, bVerbose
);
1417 ntype
= get_atomtype_ntypes(atype
);
1420 if (ir
->implicit_solvent
!= eisNO
)
1422 /* Now we have renumbered the atom types, we can check the GBSA params */
1423 check_gbsa_params(ir
,atype
);
1425 /* Check that all atoms that have charge and/or LJ-parameters also have
1426 * sensible GB-parameters
1428 check_gbsa_charges_lj(sys
,atype
);
1431 /* PELA: Copy the atomtype data to the topology atomtype list */
1432 copy_atomtype_atomtypes(atype
,&(sys
->atomtypes
));
1435 pr_symtab(debug
,0,"After renum_atype",&sys
->symtab
);
1438 fprintf(stderr
,"converting bonded parameters...\n");
1440 ntype
= get_atomtype_ntypes(atype
);
1441 convert_params(ntype
, plist
, mi
, comb
, reppow
, fudgeQQ
, sys
);
1444 pr_symtab(debug
,0,"After convert_params",&sys
->symtab
);
1446 /* set ptype to VSite for virtual sites */
1447 for(mt
=0; mt
<sys
->nmoltype
; mt
++) {
1448 set_vsites_ptype(FALSE
,&sys
->moltype
[mt
]);
1451 pr_symtab(debug
,0,"After virtual sites",&sys
->symtab
);
1453 /* Check velocity for virtual sites and shells */
1455 check_vel(sys
,state
.v
);
1461 for(i
=0; i
<sys
->nmoltype
; i
++) {
1462 check_cg_sizes(ftp2fn(efTOP
,NFILE
,fnm
),&sys
->moltype
[i
].cgs
,wi
);
1465 if (EI_DYNAMICS(ir
->eI
) && ir
->eI
!= eiBD
)
1467 check_bonds_timestep(sys
,ir
->delta_t
,wi
);
1470 check_warning_error(wi
,FARGS
);
1473 fprintf(stderr
,"initialising group options...\n");
1474 do_index(mdparin
,ftp2fn_null(efNDX
,NFILE
,fnm
),
1476 bGenVel
? state
.v
: NULL
,
1479 /* Init the temperature coupling state */
1480 init_gtc_state(&state
,ir
->opts
.ngtc
,0,ir
->opts
.nhchainlength
);
1483 fprintf(stderr
,"Checking consistency between energy and charge groups...\n");
1484 check_eg_vs_cg(sys
);
1487 pr_symtab(debug
,0,"After index",&sys
->symtab
);
1488 triple_check(mdparin
,ir
,sys
,wi
);
1489 close_symtab(&sys
->symtab
);
1491 pr_symtab(debug
,0,"After close",&sys
->symtab
);
1493 /* make exclusions between QM atoms */
1495 generate_qmexcl(sys
,ir
);
1498 if (ftp2bSet(efTRN
,NFILE
,fnm
)) {
1500 fprintf(stderr
,"getting data from old trajectory ...\n");
1501 cont_status(ftp2fn(efTRN
,NFILE
,fnm
),ftp2fn_null(efEDR
,NFILE
,fnm
),
1502 bNeedVel
,bGenVel
,fr_time
,ir
,&state
,sys
,oenv
);
1505 if (ir
->ePBC
==epbcXY
&& ir
->nwall
!=2)
1507 clear_rvec(state
.box
[ZZ
]);
1512 set_warning_line(wi
,mdparin
,-1);
1513 check_chargegroup_radii(sys
,ir
,state
.x
,wi
);
1516 if (EEL_FULL(ir
->coulombtype
)) {
1517 /* Calculate the optimal grid dimensions */
1518 copy_mat(state
.box
,box
);
1519 if (ir
->ePBC
==epbcXY
&& ir
->nwall
==2)
1520 svmul(ir
->wall_ewald_zfac
,box
[ZZ
],box
[ZZ
]);
1521 max_spacing
= calc_grid(stdout
,box
,opts
->fourierspacing
,
1522 &(ir
->nkx
),&(ir
->nky
),&(ir
->nkz
));
1523 if ((ir
->coulombtype
== eelPPPM
) && (max_spacing
> 0.1)) {
1524 set_warning_line(wi
,mdparin
,-1);
1525 warning_note(wi
,"Grid spacing larger then 0.1 while using PPPM.");
1529 if (ir
->ePull
!= epullNO
)
1530 set_pull_init(ir
,sys
,state
.x
,state
.box
,oenv
,opts
->pull_start
);
1532 /* reset_multinr(sys); */
1534 if (EEL_PME(ir
->coulombtype
)) {
1535 float ratio
= pme_load_estimate(sys
,ir
,state
.box
);
1536 fprintf(stderr
,"Estimate for the relative computational load of the PME mesh part: %.2f\n",ratio
);
1537 /* With free energy we might need to do PME both for the A and B state
1538 * charges. This will double the cost, but the optimal performance will
1539 * then probably be at a slightly larger cut-off and grid spacing.
1541 if ((ir
->efep
== efepNO
&& ratio
> 1.0/2.0) ||
1542 (ir
->efep
!= efepNO
&& ratio
> 2.0/3.0)) {
1544 "The optimal PME mesh load for parallel simulations is below 0.5\n"
1545 "and for highly parallel simulations between 0.25 and 0.33,\n"
1546 "for higher performance, increase the cut-off and the PME grid spacing");
1551 char warn_buf
[STRLEN
];
1552 double cio
= compute_io(ir
,sys
->natoms
,&sys
->groups
,F_NRE
,1);
1553 sprintf(warn_buf
,"This run will generate roughly %.0f Mb of data",cio
);
1555 set_warning_line(wi
,mdparin
,-1);
1556 warning_note(wi
,warn_buf
);
1558 printf("%s\n",warn_buf
);
1563 fprintf(stderr
,"writing run input file...\n");
1565 done_warning(wi
,FARGS
);
1567 state
.lambda
= ir
->init_lambda
;
1568 write_tpx_state(ftp2fn(efTPX
,NFILE
,fnm
),ir
,&state
,sys
);