-## Process this file with automake to produce Makefile.in
+'## Process this file with automake to produce Makefile.in
#
# Don't edit - this file is generated automatically from Makefile.am
#
ebin.h \
edsam.h \
enxio.h \
-fatal.h \
+gmx_fatal.h \
ffscanf.h \
fftgrid.h \
filenm.h \
gmx_parallel_3dfft.h \
gmx_system_xdr.h \
gmx_thread.h \
+gmx_parallel_3dfft.h \
grompp.h \
gstat.h \
index.h \
#include <stdio.h>
#include "sysstuff.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
/*
* Module to binary write and read.
#include <stdio.h>
#include "typedefs.h"
#include "main.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#define LEFT 0 /* channel to the left processor */
#define RIGHT 1 /* channel to the right processor */
#ifndef _QMMM_h
#define _QMMM_h
-static char *SRCID_QMMM_h = "$Id$";
#ifdef HAVE_IDENT
#ident "@(#) QMMM.h 1 28/2/01"
enum {
eQMmethodAM1, eQMmethodPM3, eQMmethodRHF,
- eQMmethodUHF, eQMmethodDFT, eQMmethodB3LYP, eQMmethodMP2, eQMmethodCASSCF,
+ eQMmethodUHF, eQMmethodDFT, eQMmethodB3LYP, eQMmethodMP2, eQMmethodCASSCF, eQMmethodB3LYPLAN,
eQMmethodDIRECT, eQMmethodNR
};
#include <config.h>
#endif
-#include "fatal.h"
+#include "gmx_fatal.h"
typedef enum { egcolWhite, egcolGrey, egcolBlack, egcolNR } egCol;
bool *bPerturbed;
int *typeA,*typeB;
unsigned short *ptype;
- unsigned short *cTC,*cENER,*cACC,*cFREEZE,*cVCM;
- unsigned short *cU1,*cU2,*cORF;
+ unsigned short *cTC,*cENER,*cACC,*cFREEZE,*cXTC,*cVCM;
+ unsigned short *cU1,*cU2,*cORF,*cQMMM;
/* for QMMM, atomnumber contains atomic number of the atoms */
+ int *atomnumber;
bool *bQM;
} t_mdatoms;
#include "typedefs.h"
#include "sysstuff.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#define EXP_LSB 0x00800000
#define EXP_MASK 0x7f800000
static inline void m_inv_lowerleft0(matrix src,matrix dest)
{
- double tmp = src[XX][XX]*src[YY][YY]*src[ZZ][ZZ];
- if(gmx_within_tol(tmp,0.0,100*GMX_REAL_MIN));
- gmx_fatal(FARGS,"Can not invert matrix, determinant is zero");
+ double tmp = src[XX][XX]*src[YY][YY]*src[ZZ][ZZ];
+ if (gmx_within_tol(tmp,0.0,100*GMX_REAL_MIN))
+ gmx_fatal(FARGS,"Can not invert matrix, determinant is zero");
dest[XX][XX] = 1/src[XX][XX];
dest[YY][YY] = 1/src[YY][YY];
#include "copyrite.h"
#include "sysstuff.h"
#include "txtdump.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "xtcio.h"
#include "enxio.h"
#include "smalloc.h"
#include "statutil.h"
#include "copyrite.h"
#include "confio.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "xvgr.h"
#include "gstat.h"
#include "index.h"
#include "macros.h"
#include "copyrite.h"
#include "statutil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "random.h"
#include "pdbio.h"
#include "futil.h"
#include "macros.h"
#include "copyrite.h"
#include "statutil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "random.h"
#include "strdb.h"
#include "futil.h"
#include "macros.h"
#include "copyrite.h"
#include "statutil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "random.h"
#include "pdbio.h"
#include "futil.h"
#include "macros.h"
#include "names.h"
#include "stat.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "txtdump.h"
#include "futil.h"
#include "typedefs.h"
#include "typedefs.h"
#include "statutil.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "pdbio.h"
#include "macros.h"
#include "smalloc.h"
#include "tpxio.h"
#include "statutil.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "mdatoms.h"
#include "coulomb.h"
#include "force.h"
#include "xvgr.h"
#include "pbc.h"
+
+#ifdef GMX_MPI
#include "mpi.h"
+#endif
+
#include "block_tx.h"
rvec *xptr=NULL;
top.blocks[ebCGS].multinr[i] = ncg;
}
if (PAR(cr)) {
- /* Distribute the data over processors */
- MPI_Bcast(&natoms,1,MPI_INT,root,MPI_COMM_WORLD);
/* Set some variables to zero to avoid core dumps */
ir->opts.ngtc = ir->opts.ngacc = ir->opts.ngfrz = ir->opts.ngener = 0;
+#ifdef GMX_MPI
+ /* Distribute the data over processors */
+ MPI_Bcast(&natoms,1,MPI_INT,root,MPI_COMM_WORLD);
MPI_Bcast(ir,sizeof(*ir),MPI_BYTE,root,MPI_COMM_WORLD);
MPI_Bcast(&qtot,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
+#endif
+
/* Call some dedicated communication routines, master sends n-1 times */
if (MASTER(cr)) {
for(i=1; (i<cr->nnodes); i++) {
snew(charge,natoms);
snew(x,natoms);
}
+#ifdef GMX_MPI
MPI_Bcast(charge,natoms,GMX_MPI_REAL,root,MPI_COMM_WORLD);
+#endif
}
ewaldcoeff = calc_ewaldcoeff(ir->rcoulomb,ir->ewald_rtol);
mdatoms = atoms2md(stdlog,&top.atoms,ir->opts.nFreeze,ir->eI,
ir->delta_t,0,ir->opts.tau_t,FALSE,FALSE);
snew(fr,1);
- init_forcerec(stdlog,fr,ir,&top,cr,mdatoms,nsb,box,FALSE,NULL,FALSE);
+ init_forcerec(stdlog,fr,ir,&top,cr,mdatoms,nsb,box,FALSE,NULL,NULL,FALSE);
/* First do PME based on coordinates in tpr file, send them to
* other processors if needed.
if (MASTER(cr))
fprintf(stdlog,"-----\n"
"Results based on tpr file %s\n",ftp2fn(efTPX,NFILE,fnm));
+#ifdef GMX_MPI
if (PAR(cr)) {
MPI_Bcast(x[0],natoms*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
MPI_Bcast(box[0],DIM*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
MPI_Bcast(&t,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
}
+#endif
do_my_pme(stdlog,0,bVerbose,ir,x,xbuf,f,charge,qbuf,qqbuf,box,bSort,
cr,nsb,&nrnb,&(top.atoms.excl),qtot,fr,index,NULL,
bGroups ? ir->opts.ngener : 1,mdatoms->cENER);
fp = NULL;
do {
/* Send coordinates, box and time to the other nodes */
+#ifdef GMX_MPI
if (PAR(cr)) {
MPI_Bcast(x[0],natoms*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
MPI_Bcast(box[0],DIM*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
MPI_Bcast(&t,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
}
+#endif
rm_pbc(&top.idef,nsb->natoms,box,x,x);
/* Call the PME wrapper function */
do_my_pme(stdlog,t,bVerbose,ir,x,xbuf,f,charge,qbuf,qqbuf,box,bSort,cr,
bGroups ? ir->opts.ngener : 1,mdatoms->cENER);
/* Only the master processor reads more data */
if (MASTER(cr))
- bCont = read_next_x(status,&t,natoms,x,box);
+ bCont = read_next_x(status,&t,natoms,x,box);
/* Check whether we need to continue */
+#ifdef GMX_MPI
if (PAR(cr))
- MPI_Bcast(&bCont,1,MPI_INT,root,MPI_COMM_WORLD);
+ MPI_Bcast(&bCont,1,MPI_INT,root,MPI_COMM_WORLD);
+#endif
+
} while (bCont);
/* Finish I/O, close files */
#include "typedefs.h"
#include "statutil.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "xvgr.h"
#include "pdbio.h"
#include "macros.h"
#include "typedefs.h"
#include "statutil.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "xvgr.h"
#include "pdbio.h"
#include "macros.h"
#include "smalloc.h"
#include "string2.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "atomprop.h"
#include "macros.h"
#include "index.h"
#include "network.h"
#include "block_tx.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "main.h"
#include "ns.h"
#include "macros.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "mshift.h"
#include "main.h"
#include "disre.h"
fprintf(fplog,"Step %d: bonded V and dVdl for this node\n",step);
#ifdef DEBUG
- if (g)
+ if (g && debug)
p_graph(debug,"Bondage is fun",g);
#endif
f,fshift,pbc,g,x,t1,t2,t3); /* 112 */
/* 217 TOTAL */
#ifdef DEBUG
- fprintf("idih: (%d,%d,%d,%d) cp=%g, phi=%g\n",
- ai,aj,ak,al,cos_phi,phi);
+ if (debug)
+ fprintf(debug,"idih: (%d,%d,%d,%d) cp=%g, phi=%g\n",
+ ai,aj,ak,al,cos_phi,phi);
#endif
}
#include "typedefs.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "calcgrid.h"
#define facNR 6
#include "filenm.h"
#include "pdbio.h"
#include "tpxio.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "copyrite.h"
#include "filenm.h"
#include "statutil.h"
#include "vec.h"
#include "futil.h"
#include "xvgr.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "bondf.h"
#include "copyrite.h"
#include "disre.h"
#include "vec.h"
#include "futil.h"
#include "xvgr.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "bondf.h"
#include "copyrite.h"
#include "disre.h"
#include "futil.h"
#include "string2.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "gmxfio.h"
#include "enxio.h"
#include "txtdump.h"
#include "futil.h"
#include "names.h"
-#include "fftgrid.h"
#include "writeps.h"
#include "macros.h"
#include "typedefs.h"
#include "smalloc.h"
#include "string2.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "filenm.h"
#include "futil.h"
#include "wman.h"
#include "string2.h"
#include "futil.h"
#include "network.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "statutil.h"
#include <math.h>
#include "macros.h"
#include "vec.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "gstat.h"
#include "pbc.h"
#include <ctype.h>
#include <stdio.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "smalloc.h"
#include "futil.h"
int isize[], atom_id *index[],char *grpnames[])
{
char ***gnames;
- t_block *grps;
+ t_block *grps = NULL;
int *grpnr;
snew(grpnr,ngrps);
#include "typedefs.h"
#include "smalloc.h"
#include "invblock.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
atom_id *make_invblock(const t_block *block,int nr)
{
cnt += sizeof(int);
}
-void
-F77_FUNC(xdrflong,XDRFLONG)(int *xdrid, long *lp, int *ret)
-{
- *ret = xdr_long(xdridptr[*xdrid], lp);
- cnt += sizeof(long);
-}
-
-void
F77_FUNC(xdrfshort,XDRFSHORT)(int *xdrid, short *sp, int *ret)
{
*ret = xdr_short(xdridptr[*xdrid], sp);
cnt += sizeof(char);
}
-void
-F77_FUNC(xdrfulong,XDRFULONG)(int *xdrid, unsigned long *ulp, int *ret)
-{
- *ret = xdr_u_long(xdridptr[*xdrid], (u_long *)ulp);
- cnt += sizeof(unsigned long);
-}
void
F77_FUNC(xdrfushort,XDRFUSHORT)(int *xdrid, unsigned short *usp, int *ret)
#include <limits.h>
#include <time.h>
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "network.h"
#include "main.h"
#include "macros.h"
#include "string2.h"
#include "macros.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "matio.h"
#include "statutil.h"
#include <string.h>
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "vec.h"
#include "futil.h"
#include "smalloc.h"
#include "gmxfio.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "mtxio.h"
#include "mvdata.h"
#include "network.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "symtab.h"
#include "vec.h"
#include "tgroup.h"
snew(g->acc,g->ngacc);
snew(g->nFreeze,g->ngfrz);
snew(g->egp_flags,g->ngener*g->ngener);
+
nblockrx(cr,src,g->ngtc,g->nrdf);
nblockrx(cr,src,g->ngtc,g->tau_t);
nblockrx(cr,src,g->ngtc,g->ref_t);
nblockrx(cr,src,n,g->anneal_temp[i]);
}
}
+
/* QMMM stuff, see inputrec */
blockrx(cr,src,g->ngQM);
snew(g->QMmethod,g->ngQM);
#include "mvdata.h"
#include "network.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "symtab.h"
#include "main.h"
#include "typedefs.h"
const char *eQMmethod_names[eQMmethodNR+1] = {
"AM1", "PM3", "RHF",
- "UHF", "DFT", "B3LYP", "MP2", "CASSCF",
+ "UHF", "DFT", "B3LYP", "MP2", "CASSCF","B3LYPLAN",
"DIRECT", NULL
};
#endif
#include <string.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "main.h"
#include "smalloc.h"
#include "network.h"
void gmx_left_right(int nnodes,int nodeid,int *left,int *right)
{
- *left = (nodeid - 1 + nnodes) % nnodes;
+ *left = (nnodes+nodeid-1) % nnodes;
*right = (nodeid + 1) % nnodes;
}
-
.globl nb_kernel010_ia32_3dnow
.globl _nb_kernel010_ia32_3dnow
nb_kernel010_ia32_3dnow:
-
.globl nb_kernel030_ia32_3dnow
.globl _nb_kernel030_ia32_3dnow
nb_kernel030_ia32_3dnow:
##
-
-
.globl nb_kernel101_ia32_3dnow
.globl _nb_kernel101_ia32_3dnow
nb_kernel101_ia32_3dnow:
##
-
.globl nb_kernel112_ia32_3dnow
.globl _nb_kernel112_ia32_3dnow
nb_kernel112_ia32_3dnow:
#ifndef _NB_KERNEL_IA64_SINGLE_H_
#define _NB_KERNEL_IA64_SINGLE_H_
-/*! \file nb_kernel_ia64_double.h
- * \brief ia64 double precision assembly nonbonded kernels.
+/*! \file nb_kernel_ia64_single.h
+ * \brief ia64 single precision assembly nonbonded kernels.
*
* \internal
*/
#include <stdio.h>
-#include <gmx_types.h>
-#include <gmx_neighborlist.h>
-#include <gmx_nonbonded.h>
-
-/* Include kernel headers in local directory.
- * We can only have one routine in each file due to a bug
- * in the intel assembler program...
- */
-#include "nb_kernel010_ia64_double.h"
-#include "nb_kernel010nf_ia64_double.h"
-#include "nb_kernel030_ia64_double.h"
-#include "nb_kernel030nf_ia64_double.h"
-#include "nb_kernel100_ia64_double.h"
-#include "nb_kernel100nf_ia64_double.h"
-#include "nb_kernel110_ia64_double.h"
-#include "nb_kernel110nf_ia64_double.h"
-#include "nb_kernel200_ia64_double.h"
-#include "nb_kernel200nf_ia64_double.h"
-#include "nb_kernel210_ia64_double.h"
-#include "nb_kernel210nf_ia64_double.h"
-#include "nb_kernel300_ia64_double.h"
-#include "nb_kernel300nf_ia64_double.h"
-#include "nb_kernel310_ia64_double.h"
-#include "nb_kernel310nf_ia64_double.h"
-#include "nb_kernel330_ia64_double.h"
-#include "nb_kernel330nf_ia64_double.h"
-#include "nb_kernel400_ia64_double.h"
-#include "nb_kernel400nf_ia64_double.h"
-#include "nb_kernel410_ia64_double.h"
-#include "nb_kernel410nf_ia64_double.h"
-#include "nb_kernel430_ia64_double.h"
-#include "nb_kernel430nf_ia64_double.h"
-
+#include "../nb_kerneltype.h"
+void
+nb_kernel_setup_ia64_single(FILE *log,nb_kernel_t **list);
#endif
#include "names.h"
#include "main.h"
#include "xvgr.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "physics.h"
#include "force.h"
#include "bondf.h"
tab = fr->tab14.tab;
rtab2 = sqr(fr->tab14.r);
tabscale = fr->tab14.scale;
+
krf = fr->k_rf;
crf = fr->c_rf;
itype = iatoms[i++];
ai = iatoms[i++];
aj = iatoms[i++];
- gid = GID(md->cENER[ai],md->cENER[aj],ngrp);
+ gid = GID(md->cENER[ai],md->cENER[aj],ngrp);
switch (ftype) {
case F_LJ14:
#include "random.h"
#include "bondf.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "nrama.h"
#include "rmpbc.h"
#include <string.h>
#include "sysstuff.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "names.h"
#include "macros.h"
#include "nrnb.h"
#include <string.h>
#include <ctype.h>
#include "typedefs.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "statutil.h"
#include "readinp.h"
#include "smalloc.h"
#include "pbc.h"
#include "smalloc.h"
#include "txtdump.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
/*****************************************
* PERIODIC BOUNDARY CONDITIONS
#include "txtdump.h"
#include "futil.h"
#include "names.h"
-#include "fftgrid.h"
#include "writeps.h"
#include "macros.h"
#include "xvgr.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "main.h"
btot+=bytes;
bytes/=1024;
- if ((stdlog != NULL) && (bytes != 0))
- fprintf(stdlog,"%30s:%6d kb (%7d kb) [%s, line %d, nelem %d, size %d]\n",
+ if (debug && (bytes != 0))
+ fprintf(debug,"%s:%d kb (%7d kb) [%s, line %d, nelem %d, size %d]\n",
what ? what : NN,bytes,btot/1024,
file ? file : NN,line,nelem,size);
#ifdef GMX_THREAD_PTHREAD
#include <stdio.h>
#include "typedefs.h"
#include "sysstuff.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "network.h"
#include "txtdump.h"
#include "names.h"
#include "futil.h"
#include "wman.h"
#include "tpxio.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "network.h"
#include "vec.h"
#include "string2.h"
#include "futil.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "strdb.h"
bool get_a_line(FILE *fp,char line[],int n)
#include "typedefs.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "string2.h"
#include "sysstuff.h"
#include "string2.h"
#include "typedefs.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "txtdump.h"
#include "symtab.h"
#include "typedefs.h"
#include "statutil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
/* Globals for trajectory input */
typedef struct {
-#include "fatal.h"
+#include "gmx_fatal.h"
void my_func(char *msg)
{
#include "sysstuff.h"
#include "smalloc.h"
#include "string2.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "names.h"
#include "symtab.h"
case F_HARMONIC:
case F_IDIHS:
do_harm(iparams,bRead);
+ if ((ftype == F_ANGRES || ftype == F_ANGRESZ) && bRead) {
+ /* Correct incorrect storage of parameters */
+ iparams->pdihs.phiB = iparams->pdihs.phiA;
+ iparams->pdihs.cpB = iparams->pdihs.cpA;
+ }
break;
case F_FENEBONDS:
do_real(iparams->fene.bm);
#include <string.h>
#include "sysstuff.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "txtdump.h"
#include "names.h"
#include "futil.h"
#include "widget.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
typedef struct {
Widget w,other,parent;
fprintf(out,"%s\n%s\n",GromacsVersion(),mydate(buf,255));
}
if (nldesc > 0) {
- fprintf(out,"DESCRIPTION:\n\n");
+ fprintf(out,"DESCRIPTION\n-----------\n");
print_tty_formatted(out,nldesc,desc,0);
}
if (nbug > 0) {
fprintf(out,"\n");
+ fprintf(out,"KNOWN BUGS\n----------\n");
for(i=0; i<nbug; i++) {
snew(tmp,strlen(bugs[i])+3);
strcpy(tmp,"* ");
#include <stdio.h>
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "copyrite.h"
#include "writeps.h"
#include "smalloc.h"
#ifdef HAVE_LIBXML2
#include <libxml/parser.h>
#include <libxml/tree.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "string.h"
#include "futil.h"
#include "smalloc.h"
#include "smalloc.h"
#include "vec.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#define XTC_MAGIC 1995
#include "vec.h"
#include "smalloc.h"
#include "typedefs.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "topio.h"
#include "toputil.h"
#include "convparm.h"
#include "topexcl.h"
#include "symtab.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "pgutil.h"
#include "resall.h"
#include "smalloc.h"
#include "futil.h"
#include "genalg.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "random.h"
#include "txtdump.h"
#include "vec.h"
#include "vec.h"
#include "statutil.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "physics.h"
#include "calch.h"
#include "genhydro.h"
#include "copyrite.h"
#include "sysstuff.h"
#include "txtdump.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "gmxfio.h"
#include "trnio.h"
#include "xtcio.h"
#include "copyrite.h"
#include "sysstuff.h"
#include "txtdump.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "xtcio.h"
#include "enxio.h"
#include "smalloc.h"
#include "splitter.h"
#include "sortwater.h"
#include "convparm.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "index.h"
#include "gmxfio.h"
#include "trnio.h"
static void check_pairs(int nrmols,t_molinfo mi[],int ntype,t_param *nb)
{
- int i,j,jj,k,taiA,tajA,taiB,tajB,indA,indB,bLJ;
+ int i,j,jj,k,ai,aj,taiA,tajA,taiB,tajB,bLJ;
t_params *p;
for(i=0; (i<nrmols); i++) {
p = &(mi[i].plist[F_LJ14]);
for(j=jj=0; (j<p->nr); j++) {
/* Extract atom types and hence the index in the nbf matrix */
- taiA = mi[i].atoms.atom[p->param[j].a[0]].type;
- tajA = mi[i].atoms.atom[p->param[j].a[1]].type;
- indA = ntype*taiA+tajA;
- taiB = mi[i].atoms.atom[p->param[j].a[0]].typeB;
- tajB = mi[i].atoms.atom[p->param[j].a[1]].typeB;
- indB = ntype*taiB+tajB;
+ ai = p->param[j].a[0];
+ aj = p->param[j].a[1];
+ taiA = mi[i].atoms.atom[ai].type;
+ tajA = mi[i].atoms.atom[aj].type;
+ taiB = mi[i].atoms.atom[ai].typeB;
+ tajB = mi[i].atoms.atom[aj].typeB;
bLJ = FALSE;
for(k=0; (k<MAXFORCEPARAM); k++)
- bLJ = bLJ || ((nb[indA].c[k] != 0) ||
- (nb[indB].c[k] != 0));
+ bLJ = bLJ || (((nb[taiA].c[k]*nb[tajA].c[k]) != 0) ||
+ ((nb[taiB].c[k]*nb[tajB].c[k]) != 0));
if (bLJ) {
cp_param(&(p->param[jj]),&(p->param[j]));
jj++;
}
+ else if (debug) {
+ fprintf(debug,"Removed 1-4 interaction between atoms %d and %d (within mol %s)\n",
+ ai+1,aj+1,*(mi[i].name));
+ }
}
fprintf(stderr,"Removed %d 1-4 interactions for molecule %s\n",
p->nr-jj,*(mi[i].name));
#include "pppm.h"
#include "pme.h"
#include "mdatoms.h"
+#include "rmpbc.h"
#include "repl_ex.h"
#include "qmmm.h"
#include "mpelogging.h"
#include <mpi.h>
#endif
+#ifdef USE_MPE
+#include "mpe.h"
+#include "mpelogging.h"
+#endif
+
/* The following two variables and the signal_handler function
* is used from pme.c as well
*/
* energy message in do_pmeonly => receive_lrforces) */
signal(SIGTERM,signal_handler);
signal(SIGUSR1,signal_handler);
+
+ /* Initiate PPPM if necessary */
+ if (fr->eeltype == eelPPPM)
+ init_pppm(stdlog,cr,nsb,FALSE,TRUE,state->box,getenv("GMXGHAT"),inputrec);
+ if ((fr->eeltype == eelPME) || (fr->eeltype == eelPMEUSER))
+ (void) init_pme(stdlog,cr,inputrec->nkx,inputrec->nky,inputrec->nkz,
+ inputrec->pme_order,
+ /*HOMENR(nsb),*/nsb->natoms,
+ mdatoms->bChargePerturbed,
+ inputrec->bOptFFT,inputrec->ewald_geometry);
+
+ /* Make molecules whole at start of run */
+ if (fr->ePBC != epbcNONE) {
+ do_pbc_first(stdlog,state->box,fr,graph,state->x);
+ }
/* Now do whatever the user wants us to do (how flexible...) */
switch (inputrec->eI) {
/* Initial values */
init_md(cr,inputrec,&t,&t0,&state_global->lambda,&lam0,
&mynrnb,top_global,
+
nfile,fnm,&traj,&xtc_traj,&fp_ene,&fp_dgdl,&fp_field,&mdebin,grps,
force_vir,shake_vir,mdatoms,mu_tot,&bNEMD,&bSimAnn,&vcm,nsb);
debug_gmx();
copy_rvec(state->v[ii],vcopy[ii]);
}
copy_mat(state->box,boxcopy);
- }
+ }
+
/* Write start time and temperature */
start_t=print_date_and_time(log,cr->nodeid,"Started mdrun");
}
copy_mat(boxcopy,state->box);
}
-
+
if (bVsites) {
if (graph) {
/* Following is necessary because the graph may get out of sync
shift_self(graph,state->box,state->x);
}
construct_vsites(log,state->x,&mynrnb,inputrec->delta_t,state->v,
- &top->idef,graph,cr,fr->ePBC,state->box,vsitecomm);
-
+ &top->idef,graph,cr,fr->ePBC,state->box,vsitecomm);
if (graph)
unshift_self(graph,state->box,state->x);
}
-
debug_gmx();
/* Set values for invmass etc. This routine not parallellized, but hardly
* Check comments in sim_util.c
*/
do_force(log,cr,inputrec,nsb,step,&mynrnb,top,grps,
+
state->box,state->x,f,buf,mdatoms,ener,fcd,bVerbose && !PAR(cr),
state->lambda,graph,
TRUE,bNS,FALSE,TRUE,fr,mu_tot,FALSE,t,fp_field,edyn);
#include "copyrite.h"
#include "main.h"
#include "statutil.h"
+#include "smalloc.h"
#include "futil.h"
#include "smalloc.h"
#include "edsam.h"
#include "physics.h"
#include "statutil.h"
#include "tpxio.h"
-#include "fftgrid.h"
#include "copyrite.h"
const real tol = 1e-8;
#include "vec.h"
#include "statutil.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "pdbio.h"
#include "toputil.h"
#include "h_db.h"
#include "macros.h"
#include "symtab.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "pdb2top.h"
#include "topexcl.h"
#include "topdirs.h"
#include "typedefs.h"
#include "physics.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "index.h"
#include "symtab.h"
sprintf(warn_buf,"Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
warning(NULL);
}
-
+
/* FREE ENERGY */
if (ir->efep != efepNO) {
sprintf(err_buf,"The soft-core power is %d and can only be 1 or 2",
ir->opts.ngener=atoms->grps[egcENER].nr;
nuser=str_nelem(vcm,MAXPTR,ptr1);
do_numbering(atoms,nuser,ptr1,grps,gnames,egcVCM,
- restnm,FALSE,TRUE,bVerbose);
+ restnm,forward,FALSE,TRUE,bVerbose);
/* Now we have filled the freeze struct, so we can calculate NRDF */
calc_nrdf(atoms,idef,&(ir->opts),gnames,ir->nstcomm,ir->comm_mode);
#include <string.h>
#include "typedefs.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "txtdump.h"
#include "mdrun.h"
if (re->bNPT) {
fprintf(fplog,"\nRepl p");
for(i=0; i<re->nrepl; i++)
+ {
fprintf(fplog," %5.2f",re->pres[re->ind[i]]);
- if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i-1]]))
- gmx_fatal(FARGS,"The reference pressure decreases with increasing temperature");
+ }
+
+ for(i=0; i<re->nrepl; i++)
+ {
+ if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i-1]]))
+ {
+ gmx_fatal(FARGS,"The reference pressure decreases with increasing temperature");
+ }
+ }
}
fprintf(fplog,"\nRepl ");
#include "strdb.h"
#include "futil.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "symtab.h"
#include "macros.h"
#include "resall.h"
#include "vec.h"
#include "statutil.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "pdbio.h"
#include "toputil.h"
#include "h_db.h"
#include "h_db.h"
#include "string2.h"
#include "strdb.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "ter_db.h"
#include "toputil.h"
#include "grompp.h"
#include "futil.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
typedef struct {
char *ai,*aj;
#include "names.h"
#include "string2.h"
#include "symtab.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vsite_parm.h"
#include "toputil.h"
int *nsim,
t_simsystem **sims);
-extern void generate_qmexcl(t_topology *sys,t_inputrec *ir);
+
+void
+generate_qmexcl(t_topology *sys,t_inputrec *ir);
#endif /* _topio_h */
#include "toppush.h"
#include "topdirs.h"
#include "symtab.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
void generate_nbparams(int comb,int ftype,t_params *plist,t_atomtype *atype)
#include "topdirs.h"
#include "toputil.h"
#include "symtab.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include <math.h>
/* UTILITIES */
void init_atomtype (t_atomtype *at)
{
- at->nr = 0;
- at->atom = NULL;
- at->atomname = NULL;
- at->nb = NULL;
+ at->nr = 0;
+ at->atom = NULL;
+ at->atomname = NULL;
+ at->nb = NULL;
at->bondatomtype = NULL;
- at->radius = NULL;
- at->vol = NULL;
- at->surftens = NULL;
-
+ at->radius = NULL;
+ at->vol = NULL;
+ at->surftens = NULL;
+ at->atomnumber = NULL;
}
void init_bond_atomtype (t_bond_atomtype *bat)
#include "statutil.h"
#include "sysstuff.h"
#include "txtdump.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "names.h"
#include "tpxio.h"
#include "enxio.h"
int i;
fprintf(fp,"comparing atoms\n");
+
if (a2) {
cmp_int(fp,"atoms->nr",-1,a1->nr,a2->nr);
for(i=0; (i<a1->nr); i++)
#include <math.h>
#include "index.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "string2.h"
#include "sysstuff.h"
#include "smalloc.h"
#include "physics.h"
#include "index.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "string2.h"
#include "physics.h"
force.c ghat.c init.c \
mdatom.c mdebin.c minimize.c \
ns.c nsb.c nsgrid.c \
- pme.c pppm.c \
+ pme.c pppm.c fftgrid.c \
pull.c pullinit.c \
pullio.c pullutil.c \
rf_util.c shakef.c sim_util.c \
#include "typedefs.h"
#include "network.h"
#include "vec.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "physics.h"
#include "nsb.h"
#include "main.h"
#include "macros.h"
#include "physics.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "txtdump.h"
#include "nrnb.h"
#include <stdio.h>
#include "vec.h"
#include "constr.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#ifdef DEBUG
static void check_cons(FILE *fp,char *title,real x[],int OW1,int HW2,int HW3)
#include "sysstuff.h"
#include "smalloc.h"
#include "typedefs.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "string2.h"
#include "ebin.h"
#include "main.h"
#include "smalloc.h"
#include "futil.h"
#include "fftgrid.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "physics.h"
#include "coulomb.h"
if (cr && PAR(cr))
{
#ifdef GMX_MPI
+<<<<<<< fftgrid.c
gmx_parallel_3dfft_init(&grid->mpi_fft_setup,nx,ny,nz,
cr->mpi_comm_mysim);
&(grid->pfft.local_nx),
&(grid->pfft.local_y_start_after_transpose),
&(grid->pfft.local_ny_after_transpose));
+=======
+ gmx_parallel_3dfft_init(&grid->mpi_fft_setup,nx,ny,nz,MPI_COMM_WORLD);
+
+ gmx_parallel_3dfft_limits(grid->mpi_fft_setup,
+ &(grid->pfft.local_x_start),
+ &(grid->pfft.local_nx),
+ &(grid->pfft.local_y_start_after_transpose),
+ &(grid->pfft.local_ny_after_transpose));
+>>>>>>> 1.18.2.1
#else
gmx_fatal(FARGS,"Parallel FFT supported with MPI only!");
#endif
* optimized solvent
*/
- check_solvent(fp,top,fr,&top->atoms,nsb);
-
+ check_solvent(fp,top,fr,mdatoms,nsb);
if (getenv("GMX_NO_SOLV_OPT")) {
if (fp)
#include "futil.h"
#include "vec.h"
#include "physics.h"
+<<<<<<< ghat.c
#include "coulomb.h"
-#include "fftgrid.h"
+=======
+#include "shift_util.h"
#include "xvgr.h"
+>>>>>>> 1.8.4.2
+#include "fftgrid.h"
static void calc_k(rvec lll,int ix,int iy,int iz,int nx,int ny,int nz,rvec k)
{
#include "vec.h"
#include "main.h"
#include "mvdata.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "symtab.h"
#include "txtdump.h"
#include "splittop.h"
#include "macros.h"
#include "random.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "txtdump.h"
#include "typedefs.h"
#include "update.h"
#include "ns.h"
#include "pbc.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "nrnb.h"
#include "txtdump.h"
#include "macros.h"
#include "smalloc.h"
#include "nsgrid.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "network.h"
#include "futil.h"
#include "coulomb.h"
#include "fftgrid.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "pme.h"
#include "network.h"
#include "physics.h"
#include "smalloc.h"
#include "vec.h"
#include "xvgr.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "txtdump.h"
#include "network.h"
#include "nrnb.h"
#include "pppm.h"
#include "coulomb.h"
#include "mdrun.h"
-#include "fftgrid.h"
+#include "gmx_fft.h"
#include "pme.h"
+#ifdef GMX_MPI
+#include "gmx_parallel_3dfft.h"
+#endif
+
#define llim (-1)
#define ulim (1)
#define llim2 (-3)
#include "smalloc.h"
#include "typedefs.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "rdgroup.h"
#include "symtab.h"
#include "smalloc.h"
#include "typedefs.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "rdgroup.h"
#include "symtab.h"
#include "smalloc.h"
#include "typedefs.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "rdgroup.h"
#include "symtab.h"
#include "qmmm.h"
#include <stdio.h>
#include <string.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "typedefs.h"
#include <stdlib.h>
#include "qmmm.h"
#include <stdio.h>
#include <string.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "typedefs.h"
#include <stdlib.h>
}
fprintf(out,"%s%s%s",
"%subst l9999 ",qm->devel_dir,"/l9999\n");
- fprintf(out,"%s",
- "#P ");
- fprintf(out," %s",
- eQMmethod_names[qm->QMmethod]);
- if(qm->QMmethod>=eQMmethodRHF){
- fprintf(out,"/%s",
- eQMbasis_names[qm->QMbasis]);
- if(qm->QMmethod==eQMmethodCASSCF){
- /* in case of cas, how many electrons and orbitals do we need?
- */
- fprintf(out,"(%d,%d)",
- qm->CASelectrons,qm->CASorbitals);
+ if(step){
+ fprintf(out,"%s",
+ "#T ");
+ }else{
+ fprintf(out,"%s",
+ "#P ");
+ }
+ if(qm->QMmethod==eQMmethodB3LYPLAN){
+ fprintf(out," %s",
+ "B3LYP/GEN Pseudo=Read");
+ }
+ else{
+ fprintf(out," %s",
+ eQMmethod_names[qm->QMmethod]);
+
+ if(qm->QMmethod>=eQMmethodRHF){
+ fprintf(out,"/%s",
+ eQMbasis_names[qm->QMbasis]);
+ if(qm->QMmethod==eQMmethodCASSCF){
+ /* in case of cas, how many electrons and orbitals do we need?
+ */
+ fprintf(out,"(%d,%d)",
+ qm->CASelectrons,qm->CASorbitals);
+ }
}
}
if(QMMMrec->QMMMscheme==eQMMMschemenormal){
fprintf(out," %s",
"Charge ");
}
- if (step || qm->QMmethod>=eQMmethodCASSCF){
+ if (step || qm->QMmethod==eQMmethodCASSCF){
/* fetch guess from checkpoint file, always for CASSCF */
fprintf(out,"%s"," guess=read");
}
qm->xQM[i][ZZ]/BORH2NM);
#endif
}
+
+ /* Pseudo Potential and ECP are included here if selected (MEthod suffix LAN) */
+ if(qm->QMmethod==eQMmethodB3LYPLAN){
+ fprintf(out,"\n");
+ for(i=0;i<qm->nrQMatoms;i++){
+ if(qm->atomicnumberQM[i]<21){
+ fprintf(out,"%d ",i+1);
+ }
+ }
+ fprintf(out,"\n%s\n****\n",eQMbasis_names[qm->QMbasis]);
+
+ for(i=0;i<qm->nrQMatoms;i++){
+ if(qm->atomicnumberQM[i]>21){
+ fprintf(out,"%d ",i+1);
+ }
+ }
+ fprintf(out,"\n%s\n****\n\n","lanl2dz");
+
+ for(i=0;i<qm->nrQMatoms;i++){
+ if(qm->atomicnumberQM[i]>21){
+ fprintf(out,"%d ",i+1);
+ }
+ }
+ fprintf(out,"\n%s\n","lanl2dz");
+ }
+
+
+
/* MM point charge data */
if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
fprintf(stderr,"nr mm atoms in gaussian.c = %d\n",mm->nrMMatoms);
fprintf(out,"%d F\n",i+1); /* counting from 1 */
}
}
- fprintf(out,"\n");
/* MM point charges include LJ parameters in case of QM optimization
*/
for(i=0;i<mm->nrMMatoms;i++){
}
}
fprintf(out,"\n");
+
+
fclose(out);
} /* write_gaussian_input */
#include "qmmm.h"
#include <stdio.h>
#include <string.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "typedefs.h"
#include <stdlib.h>
#include "qmmm.h"
#include <stdio.h>
#include <string.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "typedefs.h"
#include <stdlib.h>
c12au = (HARTREE2KJ*AVOGADRO*pow(BORH2NM,12));
fprintf(stderr,"there we go!\n");
+ /* fill the nucnum array in the t_mdatoms struct: */
+ for(i=0;i<md->nr;i++){
+ md->atomnumber[i] = top->atomtypes.atomnumber[md->typeA[i]];
+ }
+
/* Make a local copy of the QMMMrec */
qr = fr->qr;
#include "calcmu.h"
#include "constr.h"
#include "xvgr.h"
+
#include "mpelogging.h"
#include "domdec.h"
#endif
#include "qmmm.h"
+
#define difftime(end,start) ((double)(end)-(double)(start))
void print_time(FILE *out,time_t start,int step,t_inputrec *ir)
#include "typedefs.h"
#include "splittop.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "main.h"
#include "vsite.h"
#include "typedefs.h"
#include "names.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "xvgr.h"
#include "vec.h"
#include "xvgr.h"
#include "complex.h"
#include "copyrite.h"
-#include "fftgrid.h"
+#include "gmx_fft.h"
#include "mdrun.h"
#include "main.h"
#include "statutil.h"
+#ifdef GMX_MPI
+#include "gmx_parallel_3dfft.h"
+#endif
+
+
int main(int argc,char *argv[])
{
int mmm[] = { 8, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40,
/* We must always unshift here, also if we did not shake
* x was shifted in do_force */
+
if ((inputrec->ePBC == epbcXYZ) && (graph->nnodes > 0)) {
unshift_x(graph,state->box,state->x,xprime);
if (TRICLINIC(state->box))
#include "smalloc.h"
#include "wnblist.h"
#include "nrnb.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "macros.h"
#include "futil.h"
#include "names.h"
#include "nmol.h"
#include "manager.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#define MBFLAGS /* MB_APPLMODAL | */ MB_DONTSHOW
#include "xutil.h"
#include "copyrite.h"
#include "statutil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
/* Units are meter and second */
#include "macros.h"
#include "xutil.h"
#include "3dview.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "buttons.h"
#include "manager.h"
#include "nmol.h"
#include "sysstuff.h"
#include "macros.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "typedefs.h"
#include "string2.h"
#include "statutil.h"
#include "macros.h"
#include "xutil.h"
#include "3dview.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "buttons.h"
#include "manager.h"
#include "nmol.h"
#include "xutil.h"
#include "xdlg.h"
#include "xmb.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
/*****************************
*
* Helpful routines
#include <string.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "string2.h"
#include "sysstuff.h"
#include "smalloc.h"
#include <stdlib.h>
#include <string.h>
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "string2.h"
#include "smalloc.h"
#include "macros.h"
#include "x11.h"
#include "xdlg.h"
#include "xmb.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "gromacs.bm"
#include "stop.bm"
#include "info.bm"
#include "confio.h"
#include "copyrite.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "string2.h"
#include "vec.h"
#include "futil.h"
#include "gstat.h"
#include "names.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "string2.h"
#include "correl.h"
#include "force.h"
#include "main.h"
#include "filenm.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "mdrun.h"
#include "ns.h"
#include "txtdump.h"
#include "statutil.h"
#include "copyrite.h"
#include "pdbio.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "xvgr.h"
#include "matio.h"
#include "index.h"
#include "smalloc.h"
#include "string2.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "symtab.h"
void replace_atom(t_topology *top,int inr,char *anm,char *resnm,
#include "types/simple.h"
#include "smalloc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "sparsematrix.h"
#include "eigensolver.h"
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * GROningen Mixture of Alchemy and Childrens' Stories
+ * Green Red Orange Magenta Azure Cyan Skyblue
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
-#include "typedefs.h"
-#include "xdrf.h"
-#include "fatal.h"
-#include "smalloc.h"
+#include <gmx_ana.h>
-int xdr_real(XDR *xdrs,real *r)
-{
-#ifdef GMX_DOUBLE
- float f;
- int ret;
-
- f=*r;
- ret=xdr_float(xdrs,&f);
- *r=f;
-
- return ret;
-#else
- return xdr_float(xdrs,(float *)r);
-#endif
-}
-int xdr3drcoord(XDR *xdrs, real *fp, int *size, real *precision)
+/* This is just a wrapper binary.
+* The code that used to be in g_kinetics.c is now in gmx_kinetics.c,
+* where the old main function is called gmx_kinetics().
+*/
+int
+main(int argc, char *argv[])
{
-#ifdef GMX_DOUBLE
- float *ffp;
- float fprec;
- int i,ret,isize;
-
- isize=*size*DIM;
- if (isize <= 0)
- gmx_fatal(FARGS,"Don't know what to malloc for ffp, isize = %d",isize);
-
- snew(ffp,isize);
-
- for(i=0; (i<isize); i++)
- ffp[i]=fp[i];
- fprec=*precision;
- ret=xdr3dfcoord(xdrs,ffp,size,&fprec);
-
- *precision=fprec;
- for(i=0; (i<isize); i++)
- fp[i]=ffp[i];
-
- sfree(ffp);
- return ret;
-#else
- return xdr3dfcoord(xdrs,(float *)fp,size,(float *)precision);
-#endif
+ gmx_kinetics(argc,argv);
+ return 0;
}
-
-
-
-
-
+
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "pbc.h"
#include "copyrite.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "copyrite.h"
#include "futil.h"
#include "vec.h"
#include "index.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "xvgr.h"
#include "gstat.h"
#include "trnio.h"
#include "pbc.h"
#include "xvgr.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "statutil.h"
#include "index.h"
#include "confio.h"
#include "pdbio.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "gstat.h"
#include "macros.h"
bool bChi,bCorr,bSSHisto;
bool bDo_rt, bDo_oh, bDo_ot, bDo_jc ;
real dt=0, traj_t_ns;
-
+
atom_id isize,*index;
int ndih,nactdih,nf;
real **dih,*trans_frac,*aver_angle,*time;
#include "vec.h"
#include "index.h"
#include "pbc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "confio.h"
#include "pdbio.h"
static bool bCount=FALSE; /* just count */
static int axis = 2; /* normal to memb. default z */
static char *axtitle="Z";
- static int nslices = 10; /* nr of slices defined */
- static int ngrps = 0; /* nr. of groups */
+ static int nslices = 50; /* nr of slices defined */
+ static int ngrps = 1; /* nr. of groups */
static bool bSymmetrize=FALSE;
static bool bCenter=FALSE;
t_pargs pa[] = {
#include "vec.h"
#include "index.h"
#include "pbc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "gstat.h"
#include "pbc.h"
if (!bCutoff) {
/* write to output */
- fprintf(fp,"%8.3f ",t);
+ fprintf(fp,"%12.7f ",t);
for(g=0;(g<ngrps/2);g++) {
pbc_dx(&pbc,com[2*g],com[2*g+1],dx);
- fprintf(fp,"%10.5f %10.5f %10.5f %10.5f",
+ fprintf(fp,"%12.7f %12.7f %12.7f %12.7f",
norm(dx),dx[XX],dx[YY],dx[ZZ]);
}
fprintf(fp,"\n");
#include "copyrite.h"
#include "index.h"
#include "confio.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "physics.h"
#include "random.h"
"[TT]-angles[tt]. Both [TT]-box[tt] and [TT]-d[tt]",
"will center the system in the box.",
"[PAR]",
- "Option [TT]-bt[tt] determines the box type: [TT]tric[tt] is a",
- "triclinic box, [TT]cubic[tt] is a cubic box, [TT]dodecahedron[tt] is",
- "a rhombic dodecahedron and [TT]octahedron[tt] is a truncated octahedron.",
+ "Option [TT]-bt[tt] determines the box type: [TT]triclinic[tt] is a",
+ "triclinic box, [TT]cubic[tt] is a rectangular box with all sides equal",
+ "[TT]dodecahedron[tt] represents a rhombic dodecahedron and "
+ "[TT]octahedron[tt] is a truncated octahedron.",
"The last two are special cases of a triclinic box.",
"The length of the three box vectors of the truncated octahedron is the",
"shortest distance between two opposite hexagons.",
"[PAR]",
"Option [TT]-box[tt] requires only",
"one value for a cubic box, dodecahedron and a truncated octahedron.",
- "With [TT]-d[tt] and [TT]tric[tt] the size of the system in the x, y",
+ "[PAR]",
+ "With [TT]-d[tt] and a [TT]triclinic[tt] box the size of the system in the x, y",
"and z directions is used. With [TT]-d[tt] and [TT]cubic[tt],",
- "[TT]dodecahedron[tt] or [TT]octahedron[tt] the diameter of the system",
- "is used, which is the largest distance between two atoms.",
+ "[TT]dodecahedron[tt] or [TT]octahedron[tt] boxes, the dimensions are set",
+ "to the diameter of the system (largest distance between atoms) plus twice",
+ "the specified distance.",
"[PAR]",
"Option [TT]-angles[tt] is only meaningful with option [TT]-box[tt] and",
"a triclinic box and can not be used with option [TT]-d[tt].",
static rvec scale={1,1,1},newbox={0,0,0},newang={90,90,90};
static real rho=1000.0,rvdw=0.12;
static rvec center={0,0,0},translation={0,0,0},rotangles={0,0,0};
- static char *btype[]={ NULL, "tric", "cubic", "dodecahedron", "octahedron", NULL },*label="A";
+ static char *btype[]={ NULL, "triclinic", "cubic", "dodecahedron", "octahedron", NULL },*label="A";
static rvec visbox={0,0,0};
t_pargs pa[] = {
{ "-ndef", FALSE, etBOOL, {&bNDEF},
box[ZZ][ZZ] = sqrt(sqr(newbox[ZZ])
-box[ZZ][XX]*box[ZZ][XX]-box[ZZ][YY]*box[ZZ][YY]);
}
- if (bDist && TRICLINIC(box))
- fprintf(stderr,"WARNING: the box is triclinic, the minimum distance between the solute and the box might be less than %f\nYou can check this with g_mindist -pi\n",dist);
break;
case 'c':
case 'd':
if (check_box(box))
printf("\nWARNING: %s\n",check_box(box));
+ if (bDist && btype[0][0]=='t')
+ {
+ if(TRICLINIC(box))
+ {
+ printf("\nWARNING: Your box is triclinic with non-orthogonal axes. In this case, the\n"
+ "distance from the solute to a box surface along the corresponding normal\n"
+ "vector might be somewhat smaller than your specified value %f.\n"
+ "You can check the actual value with g_mindist -pi\n",dist);
+ }
+ else
+ {
+ printf("\nWARNING: No boxtype specified - distance condition applied in each dimension.\n"
+ "If the molecule rotates the actual distance will be smaller. You might want\n"
+ "to use a cubic box instead, or why not try a dodecahedron today?\n");
+ }
+ }
+
if (bIndex) {
fprintf(stderr,"\nSelect a group for output:\n");
get_index(&atoms,opt2fn_null("-n",NFILE,fnm),
#include "names.h"
#include "copyrite.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "enxio.h"
#define TIME_EXPLICIT 0
#include "string2.h"
#include "typedefs.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "smalloc.h"
#include "enxio.h"
#include <math.h>
#include "typedefs.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "string2.h"
#include "smalloc.h"
bORIRE = bORA || bORT || bODA || bODR || bODT;
bOTEN = opt2bSet("-oten",NFILE,fnm);
+ nset = 0;
+
snew(frame,2);
fp = open_enx(ftp2fn(efENX,NFILE,fnm),"r");
do_enxnms(fp,&nre,&enm);
#include "atomprop.h"
#include "names.h"
#include "vec.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "statutil.h"
#include "vec.h"
#include "gbutil.h"
#include "statutil.h"
#include "pbc.h"
#include "force.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "maths.h"
#include "macros.h"
#include "futil.h"
#include "macros.h"
#include "index.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
int gmx_genpr(int argc,char *argv[])
{
/*
* $Id$
*
+ *
* This source code is part of
*
* G R O M A C S
#include "tpxio.h"
#include "physics.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "index.h"
#include "smalloc.h"
#include "vec.h"
static void set_hb(t_hbdata *hb,int id,int ih, int ia,int frame,int ihb)
{
- unsigned int *ghptr;
+ unsigned int *ghptr=NULL;
if (ihb == hbHB)
ghptr = hb->hbmap[id][ia]->h[ih];
if ( !bBox || (ngrid[XX]<3) || (ngrid[YY]<3) || (ngrid[ZZ]<3) )
for(i=0; i<DIM; i++)
ngrid[i]=1;
- printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
- ngrid[XX],ngrid[YY],ngrid[ZZ],rcut);
+ else
+ printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
+ ngrid[XX],ngrid[YY],ngrid[ZZ],rcut);
snew(grid,ngrid[ZZ]);
for (z=0; z<ngrid[ZZ]; z++) {
snew((grid)[z],ngrid[YY]);
fprintf(fp,"[ %s ]",grpnames[grp]);
for (i=0; i<isize[grp]; i++) {
fprintf(fp,(i%15)?" ":"\n");
- fprintf(fp,"%4u",index[grp][i]+1);
+ fprintf(fp," %4u",index[grp][i]+1);
}
fprintf(fp,"\n");
fprintf(fp,"[ donors_hydrogens_%s ]",grpnames[grp]);
"bonds between atoms within the shell distance from the one atom are",
"considered.[PAR]"
- "It is also possible to analyse specific hydrogen bonds with",
+ /* "It is also possible to analyse specific hydrogen bonds with",
"[TT]-sel[tt]. This index file must contain a group of atom triplets",
"Donor Hydrogen Acceptor, in the following way:[PAR]",
-
+ */
"[TT]",
"[ selected ][BR]",
" 20 21 24[BR]",
{ "-merge", FALSE, etBOOL, {&bMerge},
"H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." }
};
-
+ static char *bugs[] = {
+ "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
+ };
t_filenm fnm[] = {
{ efTRX, "-f", NULL, ffREAD },
{ efTPX, NULL, NULL, ffREAD },
{ efNDX, NULL, NULL, ffOPTRD },
- { efNDX, "-sel", "select", ffOPTRD },
+ /* { efNDX, "-sel", "select", ffOPTRD },*/
{ efXVG, "-num", "hbnum", ffWRITE },
{ efXVG, "-ac", "hbac", ffOPTWR },
{ efXVG, "-dist","hbdist", ffOPTWR },
matrix box;
real t,ccut,dist,ang;
double max_nhb,aver_nhb,aver_dist;
- int h,i,j,k,l,start,end,id,ja,ogrp;
+ int h,i,j,k,l,start,end,id,ja,ogrp,nsel;
int xi,yi,zi,ai;
int xj,yj,zj,aj,xjj,yjj,zjj;
int xk,yk,zk,ak,xkk,ykk,zkk;
ppa = add_acf_pargs(&npargs,pa);
parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,NFILE,fnm,npargs,
- ppa,asize(desc),desc,0,NULL);
+ ppa,asize(desc),desc,asize(bugs),bugs);
/* process input */
bSelected = opt2bSet("-sel",NFILE,fnm);
/* analyze selected hydrogen bonds */
printf("Select group with selected atoms:\n");
get_index(&(top.atoms),opt2fn("-sel",NFILE,fnm),
- 1,isize,index,grpnames);
- if (isize[0] % 3)
+ 1,&nsel,index,grpnames);
+ if (nsel % 3)
gmx_fatal(FARGS,"Number of atoms in group '%s' not a multiple of 3\n"
"and therefore cannot contain triplets of "
"Donor-Hydrogen-Acceptor",grpnames[0]);
bTwo=FALSE;
- for(i=0; (i<isize[0]); i+=3) {
+ for(i=0; (i<nsel); i+=3) {
int dd = index[0][i];
int hh = index[0][i+1];
int aa = index[0][i+2];
add_dh (&hb->d,dd,hh,i);
add_acc(&hb->a,aa,i);
/* Should this be here ? */
- /* add_hbond(hb,dd,aa,hh,gr0,gr0,-1,FALSE,bMerge);*/
+ snew(hb->d.dptr,top.atoms.nr);
+ snew(hb->a.aptr,top.atoms.nr);
+ add_hbond(hb,dd,aa,hh,gr0,gr0,0,FALSE,bMerge,0);
}
printf("Analyzing %d selected hydrogen bonds from '%s'\n",
- hb->nrhb,grpnames[0]);
+ isize[0],grpnames[0]);
}
else {
/* analyze all hydrogen bonds: get group(s) */
if (hb->bDAnr)
count_da_grid(ngrid, grid, hb->danr[nframes]);
- /* loop over all gridcells (xi,yi,zi) */
- /* Removed confusing macro, DvdS 27/12/98 */
- for(xi=0; (xi<ngrid[XX]); xi++)
- for(yi=0; (yi<ngrid[YY]); yi++)
- for(zi=0; (zi<ngrid[ZZ]); zi++) {
+
+ if (bSelected) {
+ int ii;
+
+ for(ii=0; (ii<nsel); ii++) {
+ int dd = index[0][i];
+ int hh = index[0][i+1];
+ int aa = index[0][i+2];
+ ihb = is_hbond(hb,ii,ii,dd,aa,rcut,ccut,x,bBox,box,
+ hbox,&dist,&ang,bDA,&h,bContact);
- /* loop over donor groups gr0 (always) and gr1 (if necessary) */
- for (grp=gr0; (grp <= (bTwo?gr1:gr0)); grp++) {
- icell=&(grid[zi][yi][xi].d[grp]);
-
- if (bTwo)
- ogrp = 1-grp;
- else
- ogrp = grp;
+ if (ihb) {
+ /* add to index if not already there */
+ /* Add a hbond */
+ add_hbond(hb,dd,aa,hh,ii,ii,nframes,FALSE,bMerge,ihb);
+ }
+ }
+ }
+ else {
+ /* loop over all gridcells (xi,yi,zi) */
+ /* Removed confusing macro, DvdS 27/12/98 */
+ for(xi=0; (xi<ngrid[XX]); xi++)
+ for(yi=0; (yi<ngrid[YY]); yi++)
+ for(zi=0; (zi<ngrid[ZZ]); zi++) {
- /* loop over all hydrogen atoms from group (grp)
- * in this gridcell (icell)
- */
- for (ai=0; (ai<icell->nr); ai++) {
- i = icell->atoms[ai];
+ /* loop over donor groups gr0 (always) and gr1 (if necessary) */
+ for (grp=gr0; (grp <= (bTwo?gr1:gr0)); grp++) {
+ icell=&(grid[zi][yi][xi].d[grp]);
+
+ if (bTwo)
+ ogrp = 1-grp;
+ else
+ ogrp = grp;
- /* loop over all adjacent gridcells (xj,yj,zj) */
- /* This is a macro!!! */
- LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid,bTric) {
- jcell=&(grid[zj][yj][xj].a[ogrp]);
- /* loop over acceptor atoms from other group (ogrp)
- * in this adjacent gridcell (jcell)
- */
- for (aj=0; (aj<jcell->nr); aj++) {
- j = jcell->atoms[aj];
+ /* loop over all hydrogen atoms from group (grp)
+ * in this gridcell (icell)
+ */
+ for (ai=0; (ai<icell->nr); ai++) {
+ i = icell->atoms[ai];
+
+ /* loop over all adjacent gridcells (xj,yj,zj) */
+ /* This is a macro!!! */
+ LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid,bTric) {
+ jcell=&(grid[zj][yj][xj].a[ogrp]);
+ /* loop over acceptor atoms from other group (ogrp)
+ * in this adjacent gridcell (jcell)
+ */
+ for (aj=0; (aj<jcell->nr); aj++) {
+ j = jcell->atoms[aj];
- if ((bSelected && (j == i)) || (!bSelected)) {
/* check if this once was a h-bond */
ihb = is_hbond(hb,grp,ogrp,i,j,rcut,ccut,x,bBox,box,
hbox,&dist,&ang,bDA,&h,bContact);
-
+
if (ihb) {
/* add to index if not already there */
/* Add a hbond */
ang*=RAD2DEG;
adist[(int)( ang/abin)]++;
rdist[(int)(dist/rbin)]++;
-
+
if (!bTwo) {
int id,ia;
if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
gmx_fatal(FARGS,"Invalid donor %d",i);
if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
- gmx_fatal(FARGS,"Invalid acceptor %d",j);
+ gmx_fatal(FARGS,"Invalid acceptor %d",j);
resdist=abs(top.atoms.atom[id].resnr-
top.atoms.atom[ia].resnr);
if (resdist >= max_hx)
hb->nhx[nframes][resdist]++;
}
}
- }
- if (bInsert && bSelected) {
- /* this has been a h-bond, or we are analyzing
- selected bonds: check for inserted */
- bool ins_d, ins_a;
- real ins_d_dist, ins_d_ang, ins_a_dist, ins_a_ang;
- int ins_d_k=0,ins_a_k=0;
-
- ins_d=ins_a=FALSE;
- ins_d_dist=ins_d_ang=ins_a_dist=ins_a_ang=1e6;
-
- /* loop over gridcells adjacent to i (xk,yk,zk) */
- LOOPGRIDINNER(xk,yk,zk,xkk,ykk,zkk,xi,yi,zi,ngrid,bTric){
- kcell=&(grid[zk][yk][xk].a[grI]);
- /* loop over acceptor atoms from ins group
- in this adjacent gridcell (kcell) */
- for (ak=0; (ak<kcell->nr); ak++) {
- k=kcell->atoms[ak];
- ihb = is_hbond(hb,grp,grI,i,k,rcut,ccut,x,
- bBox,box,hbox,&dist,&ang,bDA,&h,
- bContact);
- if (ihb == hbHB) {
- if (dist < ins_d_dist) {
- ins_d=TRUE;
- ins_d_dist=dist;
- ins_d_ang =ang ;
- ins_d_k =k ;
+ if (bInsert && bSelected) {
+ /* this has been a h-bond, or we are analyzing
+ selected bonds: check for inserted */
+ bool ins_d, ins_a;
+ real ins_d_dist, ins_d_ang, ins_a_dist, ins_a_ang;
+ int ins_d_k=0,ins_a_k=0;
+
+ ins_d=ins_a=FALSE;
+ ins_d_dist=ins_d_ang=ins_a_dist=ins_a_ang=1e6;
+
+ /* loop over gridcells adjacent to i (xk,yk,zk) */
+ LOOPGRIDINNER(xk,yk,zk,xkk,ykk,zkk,xi,yi,zi,ngrid,bTric){
+ kcell=&(grid[zk][yk][xk].a[grI]);
+ /* loop over acceptor atoms from ins group
+ in this adjacent gridcell (kcell) */
+ for (ak=0; (ak<kcell->nr); ak++) {
+ k=kcell->atoms[ak];
+ ihb = is_hbond(hb,grp,grI,i,k,rcut,ccut,x,
+ bBox,box,hbox,&dist,&ang,bDA,&h,
+ bContact);
+ if (ihb == hbHB) {
+ if (dist < ins_d_dist) {
+ ins_d=TRUE;
+ ins_d_dist=dist;
+ ins_d_ang =ang ;
+ ins_d_k =k ;
+ }
}
}
}
- }
- ENDLOOPGRIDINNER;
+ ENDLOOPGRIDINNER;
/* loop over gridcells adjacent to j (xk,yk,zk) */
- LOOPGRIDINNER(xk,yk,zk,xkk,ykk,zkk,xj,yj,zj,ngrid,bTric){
- kcell=&grid[zk][yk][xk].d[grI];
- /* loop over hydrogen atoms from ins group
- in this adjacent gridcell (kcell) */
- for (ak=0; ak<kcell->nr; ak++) {
- k = kcell->atoms[ak];
- ihb = is_hbond(hb,grI,ogrp,k,j,rcut,ccut,x,
- bBox,box,hbox,&dist,&ang,bDA,&h,
- bContact);
- if (ihb == hbHB) {
- if (dist<ins_a_dist) {
- ins_a=TRUE;
- ins_a_dist=dist;
- ins_a_ang =ang ;
- ins_a_k =k ;
+ LOOPGRIDINNER(xk,yk,zk,xkk,ykk,zkk,xj,yj,zj,ngrid,bTric){
+ kcell=&grid[zk][yk][xk].d[grI];
+ /* loop over hydrogen atoms from ins group
+ in this adjacent gridcell (kcell) */
+ for (ak=0; ak<kcell->nr; ak++) {
+ k = kcell->atoms[ak];
+ ihb = is_hbond(hb,grI,ogrp,k,j,rcut,ccut,x,
+ bBox,box,hbox,&dist,&ang,bDA,&h,
+ bContact);
+ if (ihb == hbHB) {
+ if (dist<ins_a_dist) {
+ ins_a=TRUE;
+ ins_a_dist=dist;
+ ins_a_ang =ang ;
+ ins_a_k =k ;
+ }
}
}
}
- }
- ENDLOOPGRIDINNER;
+ ENDLOOPGRIDINNER;
- {
- ihb = is_hbond(hb,grI,grI,ins_d_k,ins_a_k,rcut,ccut,x,
- bBox,box,hbox,&dist,&ang,bDA,&h,
- bContact);
- if (ins_d && ins_a && ihb) {
- /* add to hbond index if not already there */
- add_hbond(hb,ins_d_k,ins_a_k,h,grI,ogrp,
- nframes,TRUE,bMerge,ihb);
-
- /* print insertion info to file */
- /*fprintf(fpins,
- "%4g: %4u:%3.3s%4d%3.3s -> "
- "%4u:%3.3s%4d%3.3s (%4.2f,%2.0f) - "
- "%4u:%3.3s%4d%3.3s (%4.2f,%2.0f)\n",t,
- a[grIA][ins_d_k]+1,
- *top.atoms.resname[top.atoms.atom[a[grIA][ins_d_k]].resnr],
- top.atoms.atom[a[grIA][ins_d_k]].resnr+1,
- *top.atoms.atomname[a[grIA][ins_d_k]],
- a[grp+grD][i]+1,
- *top.atoms.resname[top.atoms.atom[a[grp+grD][i]].resnr],
- top.atoms.atom[a[grp+grD][i]].resnr+1,
- *top.atoms.atomname[a[grp+grD][i]],
- ins_d_dist,ins_d_ang*RAD2DEG,
- a[ogrp+grA][j]+1,
- *top.atoms.resname[top.atoms.atom[a[ogrp+grA][j]].resnr],
- top.atoms.atom[a[ogrp+grA][j]].resnr+1,
- *top.atoms.atomname[a[ogrp+grA][j]],
- ins_a_dist,ins_a_ang*RAD2DEG);*/
+ {
+ ihb = is_hbond(hb,grI,grI,ins_d_k,ins_a_k,rcut,ccut,x,
+ bBox,box,hbox,&dist,&ang,bDA,&h,
+ bContact);
+ if (ins_d && ins_a && ihb) {
+ /* add to hbond index if not already there */
+ add_hbond(hb,ins_d_k,ins_a_k,h,grI,ogrp,
+ nframes,TRUE,bMerge,ihb);
+
+ /* print insertion info to file */
+ /*fprintf(fpins,
+ "%4g: %4u:%3.3s%4d%3.3s -> "
+ "%4u:%3.3s%4d%3.3s (%4.2f,%2.0f) - "
+ "%4u:%3.3s%4d%3.3s (%4.2f,%2.0f)\n",t,
+ a[grIA][ins_d_k]+1,
+ *top.atoms.resname[top.atoms.atom[a[grIA][ins_d_k]].resnr],
+ top.atoms.atom[a[grIA][ins_d_k]].resnr+1,
+ *top.atoms.atomname[a[grIA][ins_d_k]],
+ a[grp+grD][i]+1,
+ *top.atoms.resname[top.atoms.atom[a[grp+grD][i]].resnr],
+ top.atoms.atom[a[grp+grD][i]].resnr+1,
+ *top.atoms.atomname[a[grp+grD][i]],
+ ins_d_dist,ins_d_ang*RAD2DEG,
+ a[ogrp+grA][j]+1,
+ *top.atoms.resname[top.atoms.atom[a[ogrp+grA][j]].resnr],
+ top.atoms.atom[a[ogrp+grA][j]].resnr+1,
+ *top.atoms.atomname[a[ogrp+grA][j]],
+ ins_a_dist,ins_a_ang*RAD2DEG);*/
+ }
}
}
}
- }
- } /* for aj */
- }
- ENDLOOPGRIDINNER;
- } /* for ai */
- } /* for grp */
- } /* for xi,yi,zi */
+ } /* for aj */
+ }
+ ENDLOOPGRIDINNER;
+ } /* for ai */
+ } /* for grp */
+ } /* for xi,yi,zi */
+ }
analyse_donor_props(opt2fn_null("-don",NFILE,fnm),hb,nframes,t);
if (fpnhb)
do_nhb_dist(fpnhb,hb,t);
#include "confio.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "fitahx.h"
#include "futil.h"
#include "gstat.h"
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Green Red Orange Magenta Azure Cyan Skyblue
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+#include "statutil.h"
+#include "sysstuff.h"
+#include "typedefs.h"
+#include "smalloc.h"
+#include "macros.h"
+#include "gmx_fatal.h"
+#include "vec.h"
+#include "copyrite.h"
+#include "futil.h"
+#include "readinp.h"
+#include "statutil.h"
+#include "txtdump.h"
+#include "gstat.h"
+#include "xvgr.h"
+#include "physics.h"
+#include <gsl/gsl_multimin.h>
+
+enum { epAuf, epEuf, epAfu, epEfu, epNR };
+enum { eqAif, eqEif, eqAfi, eqEfi, eqAui, eqEui, eqAiu, eqEiu, eqNR };
+static char *eep[epNR] = { "Af", "Ef", "Au", "Eu" };
+static char *eeq[eqNR] = { "Aif","Eif","Afi","Efi","Aui","Eui","Aiu","Eiu" };
+
+typedef struct {
+ int nreplica; /* Number of replicas in the calculation */
+ int nframe; /* Number of time frames */
+ int nstate; /* Number of states the system can be in, e.g. F,I,U */
+ int nparams; /* Is 2, 4 or 8 */
+ bool *bMask; /* Determine whether this replica is part of the d2 comp. */
+ bool bSum;
+ int nmask; /* Number of replicas taken into account */
+ real dt; /* Timestep between frames */
+ int j0,j1; /* Range of frames used in calculating delta */
+ real **temp,**data;
+ int **state; /* State index running from 0 (F) to nstate-1 (U) */
+ real **beta,**fcalt,**icalt;
+ real *time,*sumft,*sumit,*sumfct,*sumict;
+ real *params;
+ real *d2_replica;
+} t_remd_data;
+
+char *itoa(int i)
+{
+ static char ptr[12];
+
+ sprintf(ptr,"%d",i);
+ return ptr;
+}
+
+char *epnm(int nparams,int index)
+{
+ static char buf[32],from[8],to[8];
+ int nn,ni,ii;
+
+ range_check(index,0,nparams);
+ if ((nparams == 2) || (nparams == 4))
+ return eep[index];
+ else if ((nparams > 4) && (nparams % 4 == 0))
+ return eeq[index];
+ else
+ gmx_fatal(FARGS,"Don't know how to handle %d parameters",nparams);
+
+ return NULL;
+}
+
+bool bBack(t_remd_data *d)
+{
+ return (d->nparams > 2);
+}
+
+real is_folded(t_remd_data *d,int irep,int jframe)
+{
+ if (d->state[irep][jframe] == 0)
+ return 1.0;
+ else
+ return 0.0;
+}
+
+real is_unfolded(t_remd_data *d,int irep,int jframe)
+{
+ if (d->state[irep][jframe] == d->nstate-1)
+ return 1.0;
+ else
+ return 0.0;
+}
+
+real is_intermediate(t_remd_data *d,int irep,int jframe)
+{
+ if ((d->state[irep][jframe] == 1) && (d->nstate > 2))
+ return 1.0;
+ else
+ return 0.0;
+}
+
+void integrate_dfdt(t_remd_data *d)
+{
+ int i,j;
+ double beta,ddf,ddi,df,db,fac,sumf,sumi;
+
+ d->sumfct[0] = 0;
+ d->sumict[0] = 0;
+ for(i=0; (i<d->nreplica); i++) {
+ if (d->bMask[i]) {
+ ddf = 0.5*d->dt*is_folded(d,i,0);
+ ddi = 0.5*d->dt*is_intermediate(d,i,0);
+ d->fcalt[i][0] = ddf;
+ d->icalt[i][0] = ddi;
+ d->sumfct[0] += ddf;
+ d->sumict[0] += ddi;
+ }
+ }
+ for(j=1; (j<d->nframe); j++) {
+ if (j==d->nframe-1)
+ fac = 0.5*d->dt;
+ else
+ fac = d->dt;
+ sumf = sumi = 0;
+ for(i=0; (i<d->nreplica); i++) {
+ if (d->bMask[i]) {
+ beta = d->beta[i][j];
+ if (d->nstate <= 2) {
+ df = (d->params[epAuf]*exp(-beta*d->params[epEuf])*
+ is_unfolded(d,i,j));
+ if (bBack(d)) {
+ db = (d->params[epAfu]*exp(-beta*d->params[epEfu])*
+ is_folded(d,i,j));
+ ddf = fac*(df-db);
+ }
+ else
+ ddf = fac*df;
+ d->fcalt[i][j] = d->fcalt[i][j-1] + ddf;
+ sumf += ddf;
+ }
+ else {
+ ddf = fac*((d->params[eqAif]*exp(-beta*d->params[eqEif])*
+ is_intermediate(d,i,j)) -
+ (d->params[eqAfi]*exp(-beta*d->params[eqEfi])*
+ is_folded(d,i,j)));
+ ddi = fac*((d->params[eqAui]*exp(-beta*d->params[eqEui])*
+ is_unfolded(d,i,j)) -
+ (d->params[eqAiu]*exp(-beta*d->params[eqEiu])*
+ is_intermediate(d,i,j)));
+ d->fcalt[i][j] = d->fcalt[i][j-1] + ddf;
+ d->icalt[i][j] = d->icalt[i][j-1] + ddi;
+ sumf += ddf;
+ sumi += ddi;
+ }
+ }
+ }
+ d->sumfct[j] = d->sumfct[j-1] + sumf;
+ d->sumict[j] = d->sumict[j-1] + sumi;
+ }
+}
+
+void sum_ft(t_remd_data *d)
+{
+ int i,j;
+ double fac;
+
+ for(j=0; (j<d->nframe); j++) {
+ d->sumft[j] = 0;
+ d->sumit[j] = 0;
+ if ((j == 0) || (j==d->nframe-1))
+ fac = d->dt*0.5;
+ else
+ fac = d->dt;
+ for(i=0; (i<d->nreplica); i++) {
+ if (d->bMask[i]) {
+ d->sumft[j] += fac*is_folded(d,i,j);
+ d->sumit[j] += fac*is_intermediate(d,i,j);
+ }
+ }
+ }
+}
+
+double calc_d2(t_remd_data *d)
+{
+ int i,j;
+ double dd2,d2=0,dr2,tmp;
+
+ integrate_dfdt(d);
+
+ if (d->bSum) {
+ for(j=d->j0; (j<d->j1); j++) {
+ d2 += sqr(d->sumft[j]-d->sumfct[j]);
+ if (d->nstate > 2)
+ d2 += sqr(d->sumit[j]-d->sumict[j]);
+ }
+ }
+ else {
+ for(i=0; (i<d->nreplica); i++) {
+ dr2 = 0;
+ if (d->bMask[i]) {
+ for(j=d->j0; (j<d->j1); j++) {
+ tmp = sqr(is_folded(d,i,j)-d->fcalt[i][j]);
+ d2 += tmp;
+ dr2 += tmp;
+ if (d->nstate > 2) {
+ tmp = sqr(is_intermediate(d,i,j)-d->icalt[i][j]);
+ d2 += tmp;
+ dr2 += tmp;
+ }
+ }
+ d->d2_replica[i] = dr2/(d->j1-d->j0);
+ }
+ }
+ }
+ dd2 = (d2/(d->j1-d->j0))/d->nmask;
+
+ return dd2;
+}
+
+double my_f(const gsl_vector *v,void *params)
+{
+ t_remd_data *d = (t_remd_data *) params;
+ double penalty=0;
+ int i;
+
+ for(i=0; (i<d->nparams); i++) {
+ d->params[i] = gsl_vector_get(v, i);
+ if (d->params[i] < 0)
+ penalty += 10;
+ }
+ if (penalty > 0)
+ return penalty;
+ else
+ return calc_d2(d);
+}
+
+static void optimize_remd_parameters(FILE *fp,t_remd_data *d,int maxiter,
+ real tol)
+{
+ real size,d2;
+ size_t iter = 0;
+ int status = 0;
+ int i;
+
+ const gsl_multimin_fminimizer_type *T;
+ gsl_multimin_fminimizer *s;
+
+ gsl_vector *x,*dx;
+ gsl_multimin_function my_func;
+
+ my_func.f = &my_f;
+ my_func.n = d->nparams;
+ my_func.params = (void *) d;
+
+ /* Starting point */
+ x = gsl_vector_alloc (my_func.n);
+ for(i=0; (i<my_func.n); i++)
+ gsl_vector_set (x, i, d->params[i]);
+
+ /* Step size, different for each of the parameters */
+ dx = gsl_vector_alloc (my_func.n);
+ for(i=0; (i<my_func.n); i++)
+ gsl_vector_set (dx, i, 0.1*d->params[i]);
+
+ T = gsl_multimin_fminimizer_nmsimplex;
+ s = gsl_multimin_fminimizer_alloc (T, my_func.n);
+
+ gsl_multimin_fminimizer_set (s, &my_func, x, dx);
+ gsl_vector_free (x);
+ gsl_vector_free (dx);
+
+ printf ("%5s","Iter");
+ for(i=0; (i<my_func.n); i++)
+ printf(" %8s",epnm(my_func.n,i));
+ printf (" %8s %8s\n","NM Size","Chi2");
+
+ do {
+ iter++;
+ status = gsl_multimin_fminimizer_iterate (s);
+
+ if (status != 0)
+ gmx_fatal(FARGS,"Something went wrong in the iteration in minimizer %s",
+ gsl_multimin_fminimizer_name(s));
+
+ d2 = gsl_multimin_fminimizer_minimum(s);
+ size = gsl_multimin_fminimizer_size(s);
+ status = gsl_multimin_test_size(size,tol);
+
+ if (status == GSL_SUCCESS)
+ printf ("Minimum found using %s at:\n",
+ gsl_multimin_fminimizer_name(s));
+
+ printf ("%5d", iter);
+ for(i=0; (i<my_func.n); i++)
+ printf(" %8.2e",gsl_vector_get (s->x,i));
+ printf (" %8.2e %8.2e\n",size,d2);
+ }
+ while ((status == GSL_CONTINUE) && (iter < maxiter));
+
+ gsl_multimin_fminimizer_free (s);
+}
+
+
+static void preprocess_remd(FILE *fp,t_remd_data *d,real cutoff,real tref,
+ real ucut,bool bBack,real Euf,real Efu,
+ real Ei,real t0,real t1,bool bSum)
+{
+ int i,j,ninter;
+ real dd,tau_f,tau_u;
+
+ ninter = (ucut > cutoff) ? 1 : 0;
+ if (ninter && (ucut <= cutoff))
+ gmx_fatal(FARGS,"You have requested an intermediate but the cutoff for intermediates %f is smaller than the normal cutoff(%f)",ucut,cutoff);
+
+ if (!bBack) {
+ d->nparams = 2;
+ d->nstate = 2;
+ }
+ else {
+ d->nparams = 4*(1+ninter);
+ d->nstate = 2+ninter;
+ }
+ d->bSum = bSum;
+ snew(d->beta, d->nreplica);
+ snew(d->state,d->nreplica);
+ snew(d->bMask,d->nreplica);
+ snew(d->d2_replica,d->nreplica);
+ snew(d->sumft,d->nframe);
+ snew(d->sumit,d->nframe);
+ snew(d->sumfct,d->nframe);
+ snew(d->sumict,d->nframe);
+ snew(d->params,d->nparams);
+ snew(d->fcalt,d->nreplica);
+ snew(d->icalt,d->nreplica);
+
+ if (t0 < 0)
+ d->j0 = 0;
+ else
+ for(d->j0=0; (d->j0<d->nframe) && (d->time[d->j0] < t0); d->j0++)
+ ;
+ if (t1 < 0)
+ d->j1=d->nframe;
+ else
+ for(d->j1=0; (d->j1<d->nframe) && (d->time[d->j1] < t1); d->j1++)
+ ;
+ if ((d->j1-d->j0) < d->nparams+2)
+ gmx_fatal(FARGS,"Start (%f) or end time (%f) for fitting inconsistent. Reduce t0, increase t1 or supply more data",t0,t1);
+ fprintf(fp,"Will optimize from %g to %g\n",
+ d->time[d->j0],d->time[d->j1-1]);
+ d->nmask = d->nreplica;
+ for(i=0; (i<d->nreplica); i++) {
+ snew(d->beta[i],d->nframe);
+ snew(d->state[i],d->nframe);
+ snew(d->fcalt[i],d->nframe);
+ snew(d->icalt[i],d->nframe);
+ d->bMask[i] = TRUE;
+ for(j=0; (j<d->nframe); j++) {
+ d->beta[i][j] = 1.0/(BOLTZ*d->temp[i][j]);
+ dd = d->data[i][j];
+ if (dd <= cutoff)
+ d->state[i][j] = 0;
+ else if ((ucut > cutoff) && (dd <= ucut))
+ d->state[i][j] = 1;
+ else
+ d->state[i][j] = d->nstate-1;
+ }
+ }
+ sum_ft(d);
+
+ /* Assume forward rate constant is half the total time in this
+ * simulation and backward is ten times as long */
+ tau_f = d->time[d->nframe-1];
+ tau_u = 4*tau_f;
+ d->params[epEuf] = Euf;
+ d->params[epAuf] = exp(d->params[epEuf]/(BOLTZ*tref))/tau_f;
+ if (bBack) {
+ d->params[epEfu] = Efu;
+ d->params[epAfu] = exp(d->params[epEfu]/(BOLTZ*tref))/tau_u;
+ if (ninter >0) {
+ d->params[eqEui] = Ei;
+ d->params[eqAui] = exp(d->params[eqEui]/(BOLTZ*tref))/tau_u;
+ d->params[eqEiu] = Ei;
+ d->params[eqAiu] = exp(d->params[eqEiu]/(BOLTZ*tref))/tau_u;
+ }
+ }
+ else {
+ d->params[epAfu] = 0;
+ d->params[epEfu] = 0;
+ }
+}
+
+real tau(real A,real E,real T)
+{
+ return exp(E/(BOLTZ*T))/A;
+}
+
+real folded_fraction(t_remd_data *d,real tref)
+{
+ real tauf,taub;
+
+ tauf = tau(d->params[epAuf],d->params[epEuf],tref);
+ taub = tau(d->params[epAfu],d->params[epEfu],tref);
+
+ return (taub/(tauf+taub));
+}
+
+void print_tau(FILE *gp,t_remd_data *d,real tref)
+{
+ real tauf,taub,ddd,fff,DG,DH,TDS,Tm,Tb,Te,Fb,Fe,Fm;
+ int i,np=d->nparams;
+
+ ddd = calc_d2(d);
+ fprintf(gp,"Final value for Chi2 = %12.5e (%d replicas)\n",ddd,d->nmask);
+ tauf = tau(d->params[epAuf],d->params[epEuf],tref);
+ fprintf(gp,"%s = %12.5e %s = %12.5e (kJ/mole)\n",
+ epnm(np,epAuf),d->params[epAuf],
+ epnm(np,epEuf),d->params[epEuf]);
+ if (bBack(d)) {
+ taub = tau(d->params[epAfu],d->params[epEfu],tref);
+ fprintf(gp,"%s = %12.5e %s = %12.5e (kJ/mole)\n",
+ epnm(np,epAfu),d->params[epAfu],
+ epnm(np,epEfu),d->params[epEfu]);
+ fprintf(gp,"Equilibrium properties at T = %g\n",tref);
+ fprintf(gp,"tau_f = %8.3f ns, tau_b = %8.3f ns\n",tauf/1000,taub/1000);
+ fff = taub/(tauf+taub);
+ DG = BOLTZ*tref*log(fff/(1-fff));
+ DH = d->params[epEfu]-d->params[epEuf];
+ TDS = DH-DG;
+ fprintf(gp,"Folded fraction F = %8.3f\n",fff);
+ fprintf(gp,"Unfolding energies: DG = %8.3f DH = %8.3f TDS = %8.3f\n",
+ DG,DH,TDS);
+ Tb = 260;
+ Te = 420;
+ Fb = folded_fraction(d,Tb);
+ Fe = folded_fraction(d,Te);
+ while((Te-Tb > 0.001) && (Fm != 0.5)) {
+ Tm = 0.5*(Tb+Te);
+ Fm = folded_fraction(d,Tm);
+ if (Fm > 0.5) {
+ Fb = Fm;
+ Tb = Tm;
+ }
+ else if (Fm < 0.5) {
+ Te = Tm;
+ Fe = Fm;
+ }
+ }
+ if ((Fb-0.5)*(Fe-0.5) <= 0)
+ fprintf(gp,"Melting temperature Tm = %8.3f K\n",Tm);
+ else
+ fprintf(gp,"No melting temperature detected between 260 and 420K\n");
+ if (np > 4) {
+ char *ptr;
+ fprintf(gp,"Data for intermediates at T = %g\n",tref);
+ fprintf(gp,"%8s %10s %10s %10s\n","Name","A","E","tau");
+ for(i=0; (i<np/2); i++) {
+ tauf = tau(d->params[2*i],d->params[2*i+1],tref);
+ ptr = epnm(d->nparams,2*i);
+ fprintf(gp,"%8s %10.3e %10.3e %10.3e\n",ptr+1,
+ d->params[2*i],d->params[2*i+1],tauf/1000);
+ }
+ }
+ }
+ else {
+ fprintf(gp,"Equilibrium properties at T = %g\n",tref);
+ fprintf(gp,"tau_f = %8.3f\n",tauf);
+ }
+}
+
+void dump_remd_parameters(FILE *gp,t_remd_data *d,char *fn,char *fn2,char *rfn,
+ char *efn,char *mfn,int skip,real tref)
+{
+ FILE *fp,*hp;
+ int i,j,np=d->nparams;
+ real rhs,tauf,taub,fff,DG;
+ real *params;
+ char *leg[] = { "Measured", "Fit", "Difference" };
+ char *mleg[] = { "Folded fraction","DG (kJ/mole)"};
+ char **rleg;
+ real fac[] = { 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03 };
+#define NFAC asize(fac)
+ real d2[NFAC];
+
+ integrate_dfdt(d);
+ print_tau(gp,d,tref);
+
+ if (fn) {
+ fp = xvgropen(fn,"Optimized fit to data","Time (ps)","Fraction Folded");
+ xvgr_legend(fp,asize(leg),leg);
+ for(i=0; (i<d->nframe); i++) {
+ if ((skip <= 0) || ((i % skip) == 0))
+ fprintf(fp,"%12.5e %12.5e %12.5e\n",d->time[i],
+ d->sumft[i]/d->nmask,d->sumfct[i]/d->nmask,
+ (d->sumft[i]-d->sumfct[i])/d->nmask);
+ }
+ fclose(fp);
+ }
+ if (!d->bSum && rfn) {
+ snew(rleg,d->nreplica*2);
+ for(i=0; (i<d->nreplica); i++) {
+ snew(rleg[2*i],32);
+ snew(rleg[2*i+1],32);
+ sprintf(rleg[2*i],"\\f{4}F(t) %d",i);
+ sprintf(rleg[2*i+1],"\\f{12}F \\f{4}(t) %d",i);
+ }
+ fp = xvgropen(rfn,"Optimized fit to data","Time (ps)","Fraction Folded");
+ xvgr_legend(fp,d->nreplica*2,rleg);
+ for(j=0; (j<d->nframe); j++) {
+ if ((skip <= 0) || ((j % skip) == 0)) {
+ fprintf(fp,"%12.5e",d->time[j]);
+ for(i=0; (i<d->nreplica); i++)
+ fprintf(fp," %5f %9.2e",is_folded(d,i,j),d->fcalt[i][j]);
+ fprintf(fp,"\n");
+ }
+ }
+ fclose(fp);
+ }
+
+ if (fn2 && (d->nstate > 2)) {
+ fp = xvgropen(fn2,"Optimized fit to data","Time (ps)",
+ "Fraction Intermediate");
+ xvgr_legend(fp,asize(leg),leg);
+ for(i=0; (i<d->nframe); i++) {
+ if ((skip <= 0) || ((i % skip) == 0))
+ fprintf(fp,"%12.5e %12.5e %12.5e %12.5e\n",d->time[i],
+ d->sumit[i]/d->nmask,d->sumict[i]/d->nmask,
+ (d->sumit[i]-d->sumict[i])/d->nmask);
+ }
+ fclose(fp);
+ }
+ if (mfn) {
+ if (bBack(d)) {
+ fp = xvgropen(mfn,"Melting curve","T (K)","");
+ xvgr_legend(fp,asize(mleg),mleg);
+ for(i=260; (i<=420); i++) {
+ tauf = tau(d->params[epAuf],d->params[epEuf],1.0*i);
+ taub = tau(d->params[epAfu],d->params[epEfu],1.0*i);
+ fff = taub/(tauf+taub);
+ DG = BOLTZ*i*log(fff/(1-fff));
+ fprintf(fp,"%5d %8.3f %8.3f\n",i,fff,DG);
+ }
+ fclose(fp);
+ }
+ }
+
+ if (efn) {
+ snew(params,d->nparams);
+ for(i=0; (i<d->nparams); i++)
+ params[i] = d->params[i];
+
+ hp = xvgropen(efn,"Chi2 as a function of relative parameter",
+ "Fraction","Chi2");
+ for(j=0; (j<d->nparams); j++) {
+ /* Reset all parameters to optimized values */
+ fprintf(hp,"@type xy\n");
+ for(i=0; (i<d->nparams); i++)
+ d->params[i] = params[i];
+ /* Now modify one of them */
+ for(i=0; (i<NFAC); i++) {
+ d->params[j] = fac[i]*params[j];
+ d2[i] = calc_d2(d);
+ fprintf(gp,"%s = %12g d2 = %12g\n",epnm(np,j),d->params[j],d2[i]);
+ fprintf(hp,"%12g %12g\n",fac[i],d2[i]);
+ }
+ fprintf(hp,"&\n");
+ }
+ fclose(hp);
+ for(i=0; (i<d->nparams); i++)
+ d->params[i] = params[i];
+ sfree(params);
+ }
+ if (!d->bSum) {
+ for(i=0; (i<d->nreplica); i++)
+ fprintf(gp,"Chi2[%3d] = %8.2e\n",i,d->d2_replica[i]);
+ }
+}
+
+int gmx_kinetics(int argc,char *argv[])
+{
+ static char *desc[] = {
+ "g_kinetics reads two xvg files, each one containing data for N replicas.",
+ "The first file contains the temperature of each replica at each timestep,"
+ "and the second contains real values that can be interpreted as",
+ "an indicator for folding. If the value in the file is larger than",
+ "the cutoff it is taken to be unfolded and the other way around.[PAR]",
+ "From these data an estimate of the forward and backward rate constants",
+ "for folding is made at a reference temperature. In addition,"
+ "a theoretical melting curve and free energy as a function of temperature"
+ "are printed in an xvg file.[PAR]",
+ "The user can give a max value to be regarded as intermediate",
+ "([TT]-ucut[tt]), which, when given will trigger the use of an intermediate state",
+ "in the algorithm to be defined as those structures that have",
+ "cutoff < DATA < ucut. Structures with DATA values larger than ucut will",
+ "not be regarded as potential folders. In this case 8 parameters are optimized.[PAR]",
+ "The average fraction foled is printed in an xvg file together with the fit to it.",
+ "If an intermediate is used a further file will show the build of the intermediate and the fit to that process."
+ };
+ static int nreplica = 1;
+ static real tref = 298.15;
+ static real cutoff = 0.2;
+ static real ucut = 0.0;
+ static real Euf = 10;
+ static real Efu = 30;
+ static real Ei = 10;
+ static bool bHaveT = TRUE;
+ static real t0 = -1;
+ static real t1 = -1;
+ static real tb = 0;
+ static real te = 0;
+ static real tol = 1e-3;
+ static int maxiter = 100;
+ static int skip = 0;
+ static bool bBack = TRUE;
+ static bool bSplit = TRUE;
+ static bool bSum = TRUE;
+ t_pargs pa[] = {
+ { "-time", FALSE, etBOOL, {&bHaveT},
+ "Expect a time in the input" },
+ { "-b", FALSE, etREAL, {&tb},
+ "First time to read from set" },
+ { "-e", FALSE, etREAL, {&te},
+ "Last time to read from set" },
+ { "-bfit", FALSE, etREAL, {&t0},
+ "Time to start the fit from" },
+ { "-efit", FALSE, etREAL, {&t1},
+ "Time to end the fit" },
+ { "-T", FALSE, etREAL, {&tref},
+ "Reference temperature for computing rate constants" },
+ { "-n", FALSE, etINT, {&nreplica},
+ "Read data for # replicas" },
+ { "-cut", FALSE, etREAL, {&cutoff},
+ "Cut-off (max) value for regarding a structure as folded" },
+ { "-ucut", FALSE, etREAL, {&ucut},
+ "Cut-off (max) value for regarding a structure as intermediate (if not folded)" },
+ { "-euf", FALSE, etREAL, {&Euf},
+ "Initial guess for energy of activation for folding (kJ/mole)" },
+ { "-efu", FALSE, etREAL, {&Efu},
+ "Initial guess for energy of activation for unfolding (kJ/mole)" },
+ { "-ei", FALSE, etREAL, {&Ei},
+ "Initial guess for energy of activation for intermediates (kJ/mole)" },
+ { "-maxiter", FALSE, etINT, {&maxiter},
+ "Max number of iterations" },
+ { "-back", FALSE, etBOOL, {&bBack},
+ "Take the back reaction into account" },
+ { "-tol", FALSE, etREAL,{&tol},
+ "Absolute tolerance for convergence of the Nelder and Mead simplex algorithm" },
+ { "-skip", FALSE, etINT, {&skip},
+ "Skip points in the output xvg file" },
+ { "-split", FALSE, etBOOL,{&bSplit},
+ "Estimate error by splitting the number of replicas in two and refitting" },
+ { "-sum", FALSE, etBOOL,{&bSum},
+ "Average folding before computing chi^2" }
+ };
+#define NPA asize(pa)
+
+ FILE *fp;
+ real dt_t,dt_d;
+ int nset_t,nset_d,n_t,n_d,i;
+ char *tfile,*dfile;
+ t_remd_data remd;
+
+ t_filenm fnm[] = {
+ { efXVG, "-f", "temp", ffREAD },
+ { efXVG, "-d", "data", ffREAD },
+ { efXVG, "-o", "ft_all", ffWRITE },
+ { efXVG, "-o2", "it_all", ffOPTWR },
+ { efXVG, "-o3", "ft_repl", ffOPTWR },
+ { efXVG, "-ee", "err_est", ffOPTWR },
+ { efLOG, "-g", "remd", ffWRITE },
+ { efXVG, "-m", "melt", ffWRITE }
+ };
+#define NFILE asize(fnm)
+
+ CopyRight(stderr,argv[0]);
+ parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_BE_NICE,
+ NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
+
+ if (cutoff < 0)
+ gmx_fatal(FARGS,"cutoff should be >= 0 (rather than %f)",cutoff);
+
+ tfile = opt2fn("-f",NFILE,fnm);
+ dfile = opt2fn("-d",NFILE,fnm);
+
+ fp = ffopen(opt2fn("-g",NFILE,fnm),"w");
+
+ remd.temp = read_xvg_time(tfile,bHaveT,
+ opt2parg_bSet("-b",NPA,pa),tb,
+ opt2parg_bSet("-e",NPA,pa),te,
+ nreplica,&nset_t,&n_t,&dt_t,&remd.time);
+ printf("Read %d sets of %d points in %s, dt = %g\n\n",nset_t,n_t,tfile,dt_t);
+ sfree(remd.time);
+
+ remd.data = read_xvg_time(dfile,bHaveT,
+ opt2parg_bSet("-b",NPA,pa),tb,
+ opt2parg_bSet("-e",NPA,pa),te,
+ nreplica,&nset_d,&n_d,&dt_d,&remd.time);
+ printf("Read %d sets of %d points in %s, dt = %g\n\n",nset_d,n_d,dfile,dt_d);
+
+ if ((nset_t != nset_d) || (n_t != n_d) || (dt_t != dt_d))
+ gmx_fatal(FARGS,"Files %s and %s are inconsistent. Check log file",
+ tfile,dfile);
+ remd.nreplica = nset_d;
+ remd.nframe = n_d;
+ remd.dt = dt_d;
+ preprocess_remd(fp,&remd,cutoff,tref,ucut,bBack,Euf,Efu,Ei,t0,t1,bSum);
+
+ optimize_remd_parameters(fp,&remd,maxiter,tol);
+
+ dump_remd_parameters(fp,&remd,opt2fn("-o",NFILE,fnm),
+ opt2fn_null("-o2",NFILE,fnm),
+ opt2fn_null("-o3",NFILE,fnm),
+ opt2fn_null("-ee",NFILE,fnm),
+ opt2fn("-m",NFILE,fnm),skip,tref);
+
+ if (bSplit) {
+ printf("Splitting set of replicas in two halves\n");
+ for(i=0; (i<remd.nreplica); i++)
+ remd.bMask[i] = FALSE;
+ remd.nmask = 0;
+ for(i=0; (i<remd.nreplica); i+=2) {
+ remd.bMask[i] = TRUE;
+ remd.nmask++;
+ }
+ sum_ft(&remd);
+ optimize_remd_parameters(fp,&remd,maxiter,tol);
+ dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref);
+
+ for(i=0; (i<remd.nreplica); i++)
+ remd.bMask[i] = !remd.bMask[i];
+ remd.nmask = remd.nreplica - remd.nmask;
+
+ sum_ft(&remd);
+ optimize_remd_parameters(fp,&remd,maxiter,tol);
+ dump_remd_parameters(fp,&remd,"test2.xvg",NULL,NULL,NULL,NULL,skip,tref);
+
+ for(i=0; (i<remd.nreplica); i++)
+ remd.bMask[i] = FALSE;
+ remd.nmask = 0;
+ for(i=0; (i<remd.nreplica/2); i++) {
+ remd.bMask[i] = TRUE;
+ remd.nmask++;
+ }
+ sum_ft(&remd);
+ optimize_remd_parameters(fp,&remd,maxiter,tol);
+ dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref);
+
+ for(i=0; (i<remd.nreplica); i++)
+ remd.bMask[i] = FALSE;
+ remd.nmask = 0;
+ for(i=remd.nreplica/2; (i<remd.nreplica); i++) {
+ remd.bMask[i] = TRUE;
+ remd.nmask++;
+ }
+ sum_ft(&remd);
+ optimize_remd_parameters(fp,&remd,maxiter,tol);
+ dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref);
+ }
+ fclose(fp);
+
+ view_all(NFILE, fnm);
+
+ thanx(stderr);
+
+ return 0;
+}
+
+int main(int argc, char *argv[])
+{
+ gmx_kinetics(argc,argv);
+ return 0;
+}
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "copyrite.h"
#include "futil.h"
#include "statutil.h"
#include "copyrite.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "matio.h"
#include "xvgr.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "copyrite.h"
#include "futil.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "copyrite.h"
#include "futil.h"
}
}
+
t_complex *** rc_tensor_allocation(int x, int y, int z)
{
t_complex ***t;
#include "string2.h"
#include "vec.h"
#include "index.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "princ.h"
#include "rmpbc.h"
for(i=0; (i<teller); i++) {
if ( bSplit && i>0 && abs(time[bPrev ? freq*i : i]/time_factor())<1e-5 )
fprintf(fp,"&\n");
- fprintf(fp,"%8.4f",time[bPrev ? freq*i : i]);
+ fprintf(fp,"%12.7f",time[bPrev ? freq*i : i]);
for(j=0; (j<nrms); j++) {
- fprintf(fp," %8.4f",rls[j][i]);
+ fprintf(fp," %12.7f",rls[j][i]);
if (bAv)
rlstot+=rls[j][i];
}
for(i=0; (i<teller); i++) {
if ( bSplit && i>0 && abs(time[i])<1e-5 )
fprintf(fp,"&\n");
- fprintf(fp,"%8.4f",time[i]);
+ fprintf(fp,"%12.7f",time[i]);
for(j=0; (j<nrms); j++)
- fprintf(fp," %8.4f",rlsm[j][i]);
+ fprintf(fp," %12.7f",rlsm[j][i]);
fprintf(fp,"\n");
}
}
#include "pdbio.h"
#include "tpxio.h"
#include "pbc.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "do_fit.h"
#include "princ.h"
#include "copyrite.h"
#include "index.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "xvgr.h"
#include "gstat.h"
#include "vec.h"
#include "statutil.h"
#include "copyrite.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "pbc.h"
#include "xvgr.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "copyrite.h"
#include "futil.h"
#include <math.h>
#include "confio.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "gstat.h"
#include "macros.h"
#include "confio.h"
#include "copyrite.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "gstat.h"
#include "macros.h"
#include "pbc.h"
#include "index.h"
#include "gstat.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "writeps.h"
#include "strdb.h"
#include "statutil.h"
#include "writeps.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "smalloc.h"
#include "string2.h"
#include "matio.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "vec.h"
#include "pbc.h"
#include "copyrite.h"
#include "macros.h"
#include "string2.h"
#include "futil.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
static int calc_nftype(int FTYPE,t_idef *idef)
{
#include "maths.h"
#include "string2.h"
#include "hxprops.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "smalloc.h"
#include "readev.h"
#include <stdio.h>
#include <math.h>
#include "typedefs.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "gstat.h"
real LegendreP(real x,unsigned long m)
#include "maths.h"
#include "string2.h"
#include "hxprops.h"
-#include "fatal.h"
+#include "gmx_fatal.h"
#include "futil.h"
#include "smalloc.h"
#include "readev.h"
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff