Updated x2top related files in kernel and include directories.
Fixed bug in copyrgt.
Copyrgt-ed files in contrib directory.
#include "grompp.h"
/* File argument maybe NULL, in which case the default eemprops.dat
- * is opened from the library.
+ * is opened from the library. If eemtype != -1 eemprops with eemtype
+ * equal to eemtype will be read.
*/
-extern void *read_eemprops(char *fn);
+extern void *read_eemprops(char *fn,int eemtype);
extern void write_eemprops(FILE *fp,void *eem);
-extern int eem_getnumprops(void *eem);
+extern int name2eemtype(char *name);
-extern int eem_getindex(void *eem,char *resname,char *aname,int eemtype);
+extern int eem_get_numprops(void *eem);
+
+extern int eem_get_index(void *eem,char *resname,char *aname,int eemtype);
extern real lo_get_j00(void *eem,int index,real *wj,real qH);
extern real eem_get_j00(void *eem,char *resname,char *aname,
real *wj,real qH,int eemtype);
+extern real eem_get_elem(void *eem,int index);
+
extern real eem_get_radius(void *eem,int index);
extern real eem_get_chi0(void *eem,int index);
extern void free_matrix(double **a,int n);
-extern void mat_mult(FILE *fp,int n,int m,double **x,double **y,double **z);
+extern void matrix_multiply(FILE *fp,int n,int m,double **x,double **y,double **z);
-extern void mat_inv(FILE *fp,int n,double **a);
+extern void matrix_invert(FILE *fp,int n,double **a);
#endif
#include <stdio.h>
#include "grompp.h"
-enum { eqgNone, eqgLinear, eqgYang, eqgBultinck, eqgSM, eqgNR };
+enum { eqgNone, eqgLinear, eqgYang, eqgBultinck,
+ eqgSM1, eqgSM2, eqgSM3, eqgSM4, eqgNR };
-extern void generate_charges_sm(FILE *fp,void *eem,
+extern real generate_charges_sm(FILE *fp,void *eem,
t_atoms *atoms,rvec x[],
real tol,int maxiter,void *atomprop,
- real qtotref);
+ real qtotref,int eemtype);
-extern void assign_charge_alpha(int alg,t_atoms *atoms,rvec x[],
+extern void assign_charge_alpha(int eemtype,t_atoms *atoms,rvec x[],
t_params *bonds,real tol,real fac,int maxiter,
void *atomprop,real qtotref);
AM_CPPFLAGS = -I$(top_srcdir)/include -DGMXLIBDIR=\"$(datadir)/top\"
-LDADD = ../kernel/libgmxpreprocess@LIBSUFFIX@.la ../mdlib/libmd@LIBSUFFIX@.la ../gmxlib/libgmx@LIBSUFFIX@.la
+LDADD = ../kernel/libgmxpreprocess@LIBSUFFIX@.la ../mdlib/libmd@LIBSUFFIX@.la ../gmxlib/libgmx@LIBSUFFIX@.la ../tools/libgmxana@LIBSUFFIX@.la
EXTRA_DIST = README
# issue "make <program>" and copy the binary to the correct location yourself!
# Add new entries in Makefile.am!
-EXTRA_PROGRAMS = sigeps copyrgt mkyaw \
+EXTRA_PROGRAMS = copyrgt mkyaw \
prfn options hrefify \
mkice total hexamer \
optwat addquote do_shift \
g_anavel anaf tune_dip \
- ehole pmetest gen_table
+ ehole pmetest gen_table \
+ tune_pol
# Note: you don't have to list sources for "prog" if it is the single file prog.c
pmetest_SOURCES = pmetest.c
+tune_pol_SOURCES = tune_pol.c tune_pol.h pol_xml.c poldata.c \
+ merge_xml.c
+
CLEANFILES = *~ \\\#*
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Giving Russians Opium May Alter Current Situation
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Giving Russians Opium May Alter Current Situation
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Giving Russians Opium May Alter Current Situation
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
fprintf(out,"%s %s\n",ccont,CopyrightText[i]);
for(i=0; (i<NH2); i++)
fprintf(out,"%s %s\n",ccont,head2[i]);
-
- fprintf(out,"%s %s\n",ccont,bromacs(buf,STRLEN-1));
+ bromacs(buf,STRLEN-1);
+ fprintf(out,"%s %s\n",ccont,buf);
fprintf(out,"%s\n",cend);
if (bH) {
fprintf(out,"\n");
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
extern real get_omega(real ekin,int *seed,FILE *fp,char *fn);
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
-/* $Id$ */
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.3.99_development_20071104
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2006, 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:
+ * Groningen Machine for Chemical Simulation
+ */
#include "math.h"
#include "string.h"
#include "copyrite.h"
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.3.99_development_20071104
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2006, 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:
+ * Groningen Machine for Chemical Simulation
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "tune_pol.h"
+
+static int comp_mp(const void *a,const void *b)
+{
+ t_molprop *ma = (t_molprop *)a;
+ t_molprop *mb = (t_molprop *)b;
+
+ return strcasecmp(ma->molname,mb->molname);
+}
+
+static int comp_mp_formula(const void *a,const void *b)
+{
+ int r;
+
+ t_molprop *ma = (t_molprop *)a;
+ t_molprop *mb = (t_molprop *)b;
+
+ r = strcasecmp(ma->formula,mb->formula);
+ if (r == 0)
+ return strcasecmp(ma->molname,mb->molname);
+ else
+ return r;
+}
+
+static void copy_molprop(t_molprop *dst,t_molprop *src)
+{
+ int i;
+ char *ptr;
+
+ memcpy(dst,src,sizeof(*src));
+ dst->molname = strdup(src->molname);
+ while ((ptr = strchr(dst->molname,' ')) != NULL)
+ *ptr = '-';
+ dst->formula = strdup(src->formula);
+ snew(dst->experiment,dst->nexperiment);
+ snew(dst->reference,dst->nexperiment);
+ snew(dst->pname,dst->nexperiment);
+ for(i=0; (i<dst->nexperiment); i++) {
+ dst->experiment[i] = src->experiment[i];
+ dst->reference[i] = strdup(src->reference[i]);
+ dst->pname[i] = strdup(src->pname[i]);
+ }
+}
+
+static void merge_molprop(t_molprop *dst,t_molprop *src)
+{
+ int i,nd,ns;
+
+ srenew(dst->experiment,dst->nexperiment+src->nexperiment);
+ srenew(dst->reference,dst->nexperiment+src->nexperiment);
+ srenew(dst->pname,dst->nexperiment+src->nexperiment);
+ for(i=dst->nexperiment; (i<dst->nexperiment+src->nexperiment); i++) {
+ dst->experiment[i] = src->experiment[i-dst->nexperiment];
+ dst->reference[i] = strdup(src->reference[i-dst->nexperiment]);
+ dst->pname[i] = strdup(src->pname[i-dst->nexperiment]);
+ }
+ dst->nexperiment+=src->nexperiment;
+ /* Compare src and dst for composition etc. */
+ nd = ns = 0;
+ for(i=0; (i<eatNR+eatExtra); i++) {
+ if (src->frag_comp[i] > 0)
+ ns++;
+ if (dst->frag_comp[i] > 0)
+ nd++;
+ }
+ for(i=0; (i<eelemNR); i++) {
+ if (src->elem_comp[i] > 0)
+ ns++;
+ if (dst->elem_comp[i] > 0)
+ nd++;
+ }
+ for(i=0; (i<emlNR); i++) {
+ if (src->emil_comp[i] > 0)
+ ns++;
+ if (dst->emil_comp[i] > 0)
+ nd++;
+ }
+ if (ns > 0) {
+ if (nd == 0) {
+ for(i=0; (i<eatNR+eatExtra); i++)
+ dst->frag_comp[i] = src->frag_comp[i];
+ for(i=0; (i<eelemNR); i++)
+ dst->elem_comp[i] = src->elem_comp[i];
+ for(i=0; (i<emlNR); i++)
+ dst->emil_comp[i] = src->emil_comp[i];
+ }
+ else
+ printf("Both src and dst for %s contain composition entries. Not changing anything\n",dst->molname);
+ }
+}
+
+static void clear_molprop(t_molprop *dst)
+{
+ int i;
+
+ for(i=0; (i<dst->nexperiment); i++) {
+ sfree(dst->reference[i]);
+ sfree(dst->pname[i]);
+ }
+ if (dst->nexperiment > 0) {
+ sfree(dst->experiment);
+ sfree(dst->reference);
+ sfree(dst->pname);
+ }
+ dst->nexperiment = 0;
+}
+
+static void merge_doubles(int *np,t_molprop mp[])
+{
+ int i,j,ndouble=0;
+ FILE *fp;
+
+ fp = fopen("doubles.dat","w");
+ for(i=1; (i<*np); i++) {
+ if (strcasecmp(mp[i].molname,mp[i-1].molname) == 0) {
+ if (strcasecmp(mp[i].formula,mp[i-1].formula) == 0) {
+ fprintf(fp,"%5d %s\n",ndouble+1,mp[i-1].molname);
+ merge_molprop(&(mp[i-1]),&(mp[i]));
+ for(j=i+1; (j<*np); j++) {
+ clear_molprop(&(mp[j-1]));
+ copy_molprop(&(mp[j-1]),&(mp[j]));
+ }
+ ndouble++;
+ (*np)--;
+ }
+ else {
+ printf("Molecules %s, %s have formulae %s resp. %s\n",
+ mp[i].molname,mp[i-1].molname,
+ mp[i].formula,mp[i-1].formula);
+ }
+
+ }
+ }
+ fclose(fp);
+ printf("There were %d double entries\n",ndouble);
+}
+
+static void dump_mp(int np,t_molprop mp[])
+{
+ FILE *fp;
+ int i,j,k;
+
+ fp = fopen("dump_mp.dat","w");
+
+ for(i=0; (i<np); ) {
+ for(j=i; (j<np-1) && (strcasecmp(mp[i].formula,mp[j+1].formula) == 0); j++)
+ ;
+ if (j > i) {
+ for(k=i; (k<=j); k++)
+ fprintf(fp,"%-20s %s\n",mp[k].formula,mp[k].molname);
+ fprintf(fp,"\n");
+ }
+ i=j+1;
+ }
+
+ fclose(fp);
+}
+
+t_molprop *merge_xml(int argc,char *argv[],char *outf,
+ char *sorted,int *nmolprop)
+{
+ t_molprop *mp=NULL,*mpout=NULL;
+ int i,j,np,npout=0;
+ char buf[100];
+
+ for(i=1; (i<argc); i++) {
+ np = read_molprops(argv[i],&mp,1);
+ srenew(mpout,npout+np);
+ for(j=0; (j<np); j++)
+ copy_molprop(&(mpout[npout+j]),&(mp[j]));
+ npout += np;
+ }
+
+ qsort(mpout,npout,sizeof(mp[0]),comp_mp);
+ merge_doubles(&npout,mpout);
+
+ if (outf) {
+ printf("There are %d entries to store in output file %s\n",npout,outf);
+ write_molprops(outf,npout,mpout);
+ }
+ if (sorted) {
+ qsort(mpout,npout,sizeof(mp[0]),comp_mp_formula);
+ write_molprops(sorted,npout,mpout);
+ dump_mp(npout,mpout);
+ }
+
+ *nmolprop = npout;
+
+ return mpout;
+}
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.3.99_development_20071104
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2006, 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:
+ * Groningen Machine for Chemical Simulation
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <string.h>
+#ifdef HAVE_LIBXML2
+#include <libxml/parser.h>
+#include <libxml/tree.h>
+#endif
+#include "macros.h"
+#include "tune_pol.h"
+
+extern int xmlDoValidityCheckingDefaultValue;
+
+typedef struct {
+ char *miller;
+ char *spoel;
+} t_miller2spoel;
+
+#define NM2S 4
+
+static t_miller2spoel miller2spoel[NM2S] = {
+ { "CBR", "C20" },
+ { "OPI2", "O20" },
+ { "SPI2", "S20" },
+ { "NPI2", "N20" }
+};
+
+typedef struct {
+ char *pname;
+ double pvalue,perror;
+ char *psource;
+ char *preference;
+} t_property;
+
+typedef struct {
+ char *cname;
+ int cnumber;
+} t_catom;
+
+typedef struct {
+ char *compname;
+ int ncatom;
+ t_catom *catom;
+} t_composition;
+
+typedef struct {
+ char *formula,*molname,*reference;
+ double weight;
+ int nproperty;
+ t_property *property;
+ int ncomposition;
+ t_composition *composition;
+ int ncategory;
+ char **category;
+} t_molprop2;
+
+typedef struct {
+ int npd;
+ t_molprop **pd;
+ t_molprop2 **pd2;
+} t_xmlrec;
+
+static const char *xmltypes[] = {
+ NULL,
+ "XML_ELEMENT_NODE",
+ "XML_ATTRIBUTE_NODE",
+ "XML_TEXT_NODE",
+ "XML_CDATA_SECTION_NODE",
+ "XML_ENTITY_REF_NODE",
+ "XML_ENTITY_NODE",
+ "XML_PI_NODE",
+ "XML_COMMENT_NODE",
+ "XML_DOCUMENT_NODE",
+ "XML_DOCUMENT_TYPE_NODE",
+ "XML_DOCUMENT_FRAG_NODE",
+ "XML_NOTATION_NODE",
+ "XML_HTML_DOCUMENT_NODE",
+ "XML_DTD_NODE",
+ "XML_ELEMENT_DECL",
+ "XML_ATTRIBUTE_DECL",
+ "XML_ENTITY_DECL",
+ "XML_NAMESPACE_DECL",
+ "XML_XINCLUDE_START",
+ "XML_XINCLUDE_END"
+};
+#define NXMLTYPES asize(xmltypes)
+
+enum {
+ exmlMOLECULES,
+ exmlMOLECULE, exmlFORMULA, exmlMOLNAME, exmlWEIGHT,
+ exmlCATEGORY, exmlCATNAME,
+ exmlPROPERTY, exmlPNAME, exmlPVALUE, exmlPERROR,
+ exmlPSOURCE, exmlPREFERENCE,
+ exmlCOMPOSITION, exmlCOMPNAME, exmlCATOM, exmlC_NAME, exmlC_NUMBER,
+ exmlNR
+};
+
+static const char *exml_names[exmlNR] = {
+ "molecules",
+ "molecule", "formula", "molname", "weight",
+ "category", "catname",
+ "property", "pname", "pvalue", "perror", "psource", "preference",
+ "composition", "compname", "catom", "cname", "cnumber"
+};
+
+void fatal_error(char *str,char *val)
+{
+ fprintf(stderr,"Fatal: %s - %s\n",str,val);
+ exit(1);
+}
+
+static int find_elem(char *name,int nr,const char *names[])
+{
+ int i;
+
+ for(i=0; (i<nr); i++)
+ if (strcmp(name,names[i]) == 0)
+ break;
+ if (i == nr)
+ fatal_error("Unknown element name %s",name);
+
+ return i;
+}
+
+static char *sp(int n, char buf[], int maxindent)
+{
+ int i;
+ if(n>=maxindent)
+ n=maxindent-1;
+
+ /* Don't indent more than maxindent characters */
+ for(i=0; (i<n); i++)
+ buf[i] = ' ';
+ buf[i] = '\0';
+
+ return buf;
+}
+
+static void set_miller(t_molprop *mp,char *milnm,int number,int nh)
+{
+ int i;
+
+ mp->emil_comp[emlH] += nh*number;
+ for(i=0; (i<emlNR); i++) {
+ if (strcmp(miller[i].name,milnm) == 0) {
+ mp->emil_comp[i] += number;
+ return;
+ }
+ }
+ if (i == emlNR)
+ fatal_error(">No such Miller type",milnm);
+}
+
+static void set_bosque(t_molprop *mp,char *bosnm,int number,int nh)
+{
+ int i;
+
+ mp->elem_comp[eelemH] += nh*number;
+ for(i=0; (i<eelemNR); i++) {
+ if (strcmp(bosque[i].name,bosnm) == 0) {
+ mp->elem_comp[i] += number;
+ return;
+ }
+ }
+ if (i == eelemNR)
+ fatal_error(">No such Bosque type",bosnm);
+}
+
+static void copy_molprop(t_molprop *dst,t_molprop2 *src,int update_bm)
+{
+ int i,j,k,n,nh,imiller;
+ char *ptr;
+
+ if (src->formula)
+ dst->formula = strdup(src->formula);
+ if (src->molname)
+ dst->molname = strdup(src->molname);
+ dst->weight = src->weight;
+
+ for(i=0; (i<src->ncategory); i++) {
+ for(j=0; (j<ecNR); j++) {
+ if (strcasecmp(src->category[i],ec_name[j]) == 0) {
+ dst->category[j] = 1;
+ }
+ }
+ }
+
+ for(i=0;(i<src->nproperty);i++){
+ if (strcasecmp(src->property[i].psource,"experiment") == 0) {
+ dst->nexperiment++;
+ srenew(dst->experiment,dst->nexperiment);
+ srenew(dst->reference,dst->nexperiment);
+ srenew(dst->pname,dst->nexperiment);
+ dst->experiment[dst->nexperiment-1] = src->property[i].pvalue;
+ dst->reference[dst->nexperiment-1] = strdup(src->property[i].preference);
+ dst->pname[dst->nexperiment-1] = strdup(src->property[i].pname);
+ }
+ else {
+ for(j=0; (j<eqmNR); j++)
+ if (strcasecmp(lbasis[j],src->property[i].psource) == 0)
+ break;
+ if (j == eqmNR)
+ fatal_error("No such method known",src->property[i].psource);
+ dst->qm[j] = src->property[i].pvalue;
+ }
+ }
+
+ for(i=0; (i<src->ncomposition); i++) {
+ if (strcasecmp(src->composition[i].compname,"spoel") == 0) {
+ for(k=0; (k<src->composition[i].ncatom); k++) {
+ ptr = src->composition[i].catom[k].cname;
+ n = src->composition[i].catom[k].cnumber;
+
+ for(imiller=0; (imiller<NM2S); imiller++) {
+ if (strcasecmp(miller2spoel[imiller].miller,ptr) == 0) {
+ ptr = miller2spoel[imiller].spoel;
+ if (update_bm)
+ set_miller(dst,miller2spoel[imiller].miller,n,0);
+ printf("Converting %s to %s\n",miller2spoel[imiller].miller,
+ miller2spoel[imiller].spoel);
+ break;
+ }
+ }
+
+ for(j=0; (j<eatNR+eatExtra); j++)
+ if (strcasecmp(spoel[j].name,ptr) == 0) {
+ dst->frag_comp[j] += n;
+ if (update_bm) {
+ nh = spoel[j].nh;
+ if (imiller == NM2S)
+ set_miller(dst,spoel[j].miller,n,nh);
+ set_bosque(dst,spoel[j].bosque,n,nh);
+ }
+ break;
+ }
+ }
+ }
+ else if ((strcasecmp(src->composition[i].compname,"bosque") == 0) && !update_bm) {
+ for(k=0; (k<src->composition[i].ncatom); k++) {
+ for(j=0; (j<eelemNR); j++)
+ if (strcasecmp(bosque[j].name,src->composition[i].catom[k].cname) == 0) {
+ dst->elem_comp[j] += src->composition[i].catom[k].cnumber;
+ break;
+ }
+ }
+ }
+ else if ((strcasecmp(src->composition[i].compname,"miller") == 0) && !update_bm) {
+ for(k=0; (k<src->composition[i].ncatom); k++) {
+ for(j=0; (j<emlNR); j++)
+ if (strcasecmp(miller[j].name,src->composition[i].catom[k].cname) == 0) {
+ dst->emil_comp[j] += src->composition[i].catom[k].cnumber;
+ break;
+ }
+ }
+ }
+ }
+}
+
+static void process_attr(FILE *fp,xmlAttrPtr attr,int elem,
+ int indent,t_xmlrec *xml)
+ {
+ char *attrname,*attrval;
+ char buf[100];
+ t_molprop2 *xptr=NULL;
+ t_catom *ca;
+
+ if (xml->npd > 0)
+ xptr = xml->pd2[xml->npd-1];
+
+ while (attr != NULL) {
+ attrname = (char *)attr->name;
+ attrval = (char *)attr->children->content;
+
+#define atest(s) ((strcasecmp(attrname,s) == 0) && (attrval != NULL))
+ switch (elem) {
+ case exmlMOLECULE:
+ if (atest("formula"))
+ xptr->formula = strdup(attrval);
+ else if (atest("molname"))
+ xptr->molname = strdup(attrval);
+ else if (atest("weight"))
+ xptr->weight = atof(attrval);
+ else
+ fatal_error("Unknown attribute for molecule",attrname);
+ break;
+
+ case exmlCATEGORY:
+ if (atest("catname"))
+ xptr->category[xptr->ncategory-1] = strdup(attrval);
+ else
+ fatal_error("Unknown attribute for category",attrname);
+ break;
+
+ case exmlPROPERTY: {
+ t_property *xp = &(xptr->property[xptr->nproperty-1]);
+ if (atest("pname"))
+ xp->pname = strdup(attrval);
+ else if (atest("psource"))
+ xp->psource = strdup(attrval);
+ else if (atest("preference"))
+ xp->preference = strdup(attrval);
+ else if (atest("pvalue"))
+ xp->pvalue = atof(attrval);
+ else if (atest("perror"))
+ xp->perror = atof(attrval);
+ else
+ fatal_error("Unknown attribute for property",attrname);
+ break;
+ }
+ case exmlCOMPOSITION: {
+ if (xptr->ncomposition <= 0)
+ fatal_error("Composition error","");
+ t_composition *co = &(xptr->composition[xptr->ncomposition-1]);
+ if (atest("compname"))
+ co->compname = strdup(attrval);
+ else
+ fatal_error("Unknown attribute for catom",attrname);
+ break;
+ }
+ case exmlCATOM:
+ if (xptr->ncomposition <= 0)
+ fatal_error("Composition error","");
+ ca = &(xptr->composition[xptr->ncomposition-1].catom[xptr->composition[xptr->ncomposition-1].ncatom-1]);
+ if (atest("cname"))
+ ca->cname = strdup(attrval);
+ else if (atest("cnumber"))
+ ca->cnumber = atoi(attrval);
+ else
+ fatal_error("Unknown attribute for catom",attrname);
+ break;
+
+ default:
+ break;
+ }
+ if (fp)
+ fprintf(fp,"%sProperty: '%s' Value: '%s'\n",sp(indent,buf,99),
+ attrname,attrval);
+ attr = attr->next;
+#undef atest
+ }
+}
+
+static void process_element(FILE *fp,xmlNodePtr tree,int indent,t_xmlrec *xml)
+{
+ int elem;
+ char buf[100];
+ t_molprop2 *xptr=NULL;
+
+ if (xml->npd > 0)
+ xptr = xml->pd2[xml->npd-1];
+
+ elem = find_elem((char *)tree->name,exmlNR,exml_names);
+ if (fp)
+ fprintf(fp,"%sElement node name %s\n",sp(indent,buf,99),
+ (char *)tree->name);
+ switch (elem) {
+ case exmlMOLECULES:
+ break;
+ case exmlMOLECULE:
+ xml->npd++;
+ srenew(xml->pd2,xml->npd);
+ snew((xml->pd2[xml->npd-1]),1);
+ break;
+ case exmlCATEGORY:
+ xptr->ncategory++;
+ srenew(xptr->category,xptr->ncategory);
+ nullify(xptr->category[xptr->ncategory-1]);
+ break;
+ case exmlPROPERTY:
+ xptr->nproperty++;
+ srenew(xptr->property,xptr->nproperty);
+ nullify(xptr->property[xptr->nproperty-1]);
+ break;
+ case exmlCOMPOSITION:
+ xptr->ncomposition++;
+ srenew(xptr->composition,xptr->ncomposition);
+ nullify(xptr->composition[xptr->ncomposition-1]);
+ break;
+ case exmlCATOM: {
+ t_composition *co = &(xptr->composition[xptr->ncomposition-1]);
+
+ co->ncatom++;
+ srenew(co->catom,co->ncatom);
+ nullify(co->catom[co->ncatom-1]);
+ break;
+ }
+ default:
+ break;
+ }
+ process_attr(fp,tree->properties,elem,indent+2,xml);
+}
+
+static void process_tree(FILE *fp,xmlNodePtr tree,int indent,t_xmlrec *xml)
+{
+ char buf[100];
+
+ while (tree != NULL) {
+ if (fp) {
+ if ((tree->type > 0) && (tree->type < NXMLTYPES))
+ fprintf(fp,"Node type %s encountered with name %s\n",
+ xmltypes[tree->type],(char *)tree->name);
+ else
+ fprintf(fp,"Node type %d encountered\n",tree->type);
+ }
+
+ switch (tree->type) {
+ case XML_ELEMENT_NODE:
+ process_element(fp,tree,indent+2,xml);
+
+ if (tree->children)
+ process_tree(fp,tree->children,indent+2,xml);
+ break;
+ default:
+ break;
+ }
+ tree = tree->next;
+ }
+}
+
+static t_xmlrec *init_xml()
+{
+ t_xmlrec *xml;
+
+ snew(xml,1);
+
+ return xml;
+}
+
+static void done_xml(t_xmlrec **xml)
+{
+ free((*xml)->pd);
+ free((*xml)->pd2);
+ free(*xml);
+ *xml = NULL;
+}
+
+int read_molprops(char *fn,t_molprop **pd,int update_bm)
+{
+ xmlDocPtr doc;
+ t_xmlrec *xml;
+ int i,npd;
+
+ xmlDoValidityCheckingDefaultValue = 0;
+ if ((doc = xmlParseFile(fn)) == NULL) {
+ fprintf(stderr,"Reading XML file %s. Run a syntax checker such as nsgmls.",
+ fn);
+ exit(1);
+ }
+
+ xml = init_xml();
+ process_tree(NULL,doc->children,0,xml);
+
+ /* Copy xml to pd */
+ npd = xml->npd;
+ snew(*pd,npd);
+ for(i=0; (i<npd); i++)
+ copy_molprop(&((*pd)[i]),xml->pd2[i],update_bm);
+
+ xmlFreeDoc(doc);
+ done_xml(&xml);
+
+ return npd;
+}
+
+static void add_xml_int(xmlNodePtr ptr,xmlChar *name,int val)
+{
+ xmlChar buf[32];
+
+ sprintf((char *)buf,"%d",val);
+ if (xmlSetProp(ptr,name,buf) == 0)
+ fatal_error("Setting",(char *)name);
+}
+
+static void add_xml_double(xmlNodePtr ptr,xmlChar *name,double val)
+{
+ xmlChar buf[32];
+
+ sprintf((char *)buf,"%g",val);
+ if (xmlSetProp(ptr,name,buf) == 0)
+ fatal_error("Setting",(char *)name);
+}
+
+static void add_xml_char(xmlNodePtr ptr,xmlChar *name,char *val)
+{
+ if (xmlSetProp(ptr,name,(xmlChar *)val) == 0)
+ fatal_error("Setting",(char *)name);
+}
+
+static xmlNodePtr add_xml_child(xmlNodePtr parent,int type)
+{
+ xmlNodePtr child;
+
+ if ((child = xmlNewChild(parent,NULL,(xmlChar *)exml_names[type],NULL)) == NULL)
+ fatal_error("Creating element",(char *)exml_names[type]);
+
+ return child;
+}
+
+static xmlNodePtr add_xml_comment(xmlDocPtr doc,
+ xmlNodePtr prev,xmlChar *comment)
+{
+ xmlNodePtr comm,ptr;
+
+ if ((comm = xmlNewComment(comment)) == NULL)
+ fatal_error("Creating doc comment element","");
+ ptr = prev;
+ while (ptr->next != NULL)
+ ptr=ptr->next;
+ ptr->next = comm;
+ comm->prev = ptr;
+ comm->doc = doc;
+
+ return comm;
+}
+
+static void add_xml_molprop(xmlNodePtr parent,t_molprop *mp)
+{
+ xmlNodePtr ptr,child,comp;
+ int i;
+
+ ptr = add_xml_child(parent,exmlMOLECULE);
+ add_xml_char(ptr,(xmlChar *)exml_names[exmlMOLNAME],mp->molname);
+ add_xml_char(ptr,(xmlChar *)exml_names[exmlFORMULA],mp->formula);
+ add_xml_double(ptr,(xmlChar *)exml_names[exmlWEIGHT],mp->weight);
+
+ for(i=0; (i<ecNR); i++)
+ if (mp->category[i] > 0) {
+ child = add_xml_child(ptr,exmlCATEGORY);
+ add_xml_char(child,(xmlChar *)exml_names[exmlCATNAME],ec_name[i]);
+ }
+
+ for(i=0; (i<mp->nexperiment); i++) {
+ child = add_xml_child(ptr,exmlPROPERTY);
+ add_xml_char(child,(xmlChar *)exml_names[exmlPNAME],mp->pname[i]);
+ add_xml_double(child,(xmlChar *)exml_names[exmlPVALUE],mp->experiment[i]);
+ add_xml_char(child,(xmlChar *)exml_names[exmlPSOURCE],"Experiment");
+ add_xml_char(child,(xmlChar *)exml_names[exmlPREFERENCE],mp->reference[i]);
+ }
+ for(i=0; (i<eqmNR); i++) {
+ if (mp->qm[i] > 0) {
+ child = add_xml_child(ptr,exmlPROPERTY);
+ add_xml_char(child,(xmlChar *)exml_names[exmlPNAME],"Polarizability");
+ add_xml_double(child,(xmlChar *)exml_names[exmlPVALUE],mp->qm[i]);
+ add_xml_char(child,(xmlChar *)exml_names[exmlPSOURCE],lbasis[i]);
+ add_xml_char(child,(xmlChar *)exml_names[exmlPREFERENCE],"Spoel2007a");
+ }
+ }
+ child = add_xml_child(ptr,exmlCOMPOSITION);
+ add_xml_char(child,(xmlChar *)exml_names[exmlCOMPNAME],"Miller");
+ for(i=0; (i<emlNR); i++)
+ if (mp->emil_comp[i] > 0) {
+ comp = add_xml_child(child,exmlCATOM);
+ add_xml_char(comp,(xmlChar *)exml_names[exmlC_NAME],miller[i].name);
+ add_xml_int(comp,(xmlChar *)exml_names[exmlC_NUMBER],mp->emil_comp[i]);
+ }
+ child = add_xml_child(ptr,exmlCOMPOSITION);
+ add_xml_char(child,(xmlChar *)exml_names[exmlCOMPNAME],"Bosque");
+ for(i=0; (i<eelemNR); i++)
+ if (mp->elem_comp[i] > 0) {
+ comp = add_xml_child(child,exmlCATOM);
+ add_xml_char(comp,(xmlChar *)exml_names[exmlC_NAME],bosque[i].name);
+ add_xml_int(comp,(xmlChar *)exml_names[exmlC_NUMBER],mp->elem_comp[i]);
+ }
+ child = add_xml_child(ptr,exmlCOMPOSITION);
+ add_xml_char(child,(xmlChar *)exml_names[exmlCOMPNAME],"Spoel");
+ for(i=0; (i<eatNR+eatExtra); i++)
+ if (mp->frag_comp[i] > 0) {
+ comp = add_xml_child(child,exmlCATOM);
+ add_xml_char(comp,(xmlChar *)exml_names[exmlC_NAME],spoel[i].name);
+ add_xml_int(comp,(xmlChar *)exml_names[exmlC_NUMBER],mp->frag_comp[i]);
+ }
+}
+
+void write_molprops(char *fn,int npd,t_molprop pd[])
+{
+ xmlDocPtr doc;
+ xmlDtdPtr dtd;
+ xmlNodePtr myroot;
+ int i,nmt;
+ xmlChar *libdtdname,*dtdname,*gmx;
+
+ gmx = (xmlChar *) "molecules";
+ dtdname = (xmlChar *) "molprops.dtd";
+ libdtdname = dtdname;
+
+ if ((doc = xmlNewDoc((xmlChar *)"1.0")) == NULL)
+ fatal_error("Creating XML document","");
+
+ if ((dtd = xmlCreateIntSubset(doc,dtdname,libdtdname,dtdname)) == NULL)
+ fatal_error("Creating XML DTD","");
+
+ if ((myroot = xmlNewDocNode(doc,NULL,gmx,NULL)) == NULL)
+ fatal_error("Creating root element","");
+ dtd->next = myroot;
+ myroot->prev = (xmlNodePtr) dtd;
+
+ /* Add molecule definitions */
+ for(i=0; (i<npd); i++)
+ add_xml_molprop(myroot,&(pd[i]));
+
+ xmlSetDocCompressMode(doc,0);
+ xmlIndentTreeOutput = 1;
+ if (xmlSaveFormatFileEnc(fn,doc,"ISO-8859-1",2) == 0)
+ fatal_error("Saving file",fn);
+ xmlFreeDoc(doc);
+}
+
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.3.99_development_20071104
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2006, 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:
+ * Groningen Machine for Chemical Simulation
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "string2.h"
+#include "tune_pol.h"
+
+/* Noble gases and cations from G. D. Mahan. Modied Sternheimer
+ equation for polarizability. Phys. Rev. A, 22: 1780-1785, 1980. */
+
+/* Anions from Christof H{\"a}ttig and Bernd Artur Hess. {TDMP2}
+ calculation of dynamic multipole polarizabilities and dispersion
+ coefficients for the halogen anions F-, Cl-, Br- and
+ I-. J. Chem. Phys., 108:3863-3870, 1998. */
+
+t_spoel spoel[eatNR+eatExtra] = {
+ { "C33", "CH3", "CTE", "C", 3, 4, 3, 0, 0, 0.15 },
+ { "C32", "CH2", "CTE", "C", 2, 4, 3, 0, 0, 0.15 },
+ { "C31", "CH", "CTE", "C", 1, 4, 3, 0, 0, 0.15 },
+ { "C30", "C", "CTE", "C", 0, 4, 3, 0, 0, 0.15 },
+ { "C22", "CH2", "CTR", "C", 2, 3, 2, 0, 0, 0.14 },
+ { "C21", "CH", "CTR", "C", 1, 3, 2, 0, 0, 0.14 },
+ { "C20", "C", "CTR", "C", 0, 3, 2, 0, 0, 0.14 },
+ { "C11", "CH", "CDI", "C", 1, 2, 1, 0, 0, 0.13 },
+ { "C10", "C", "CDI", "C", 0, 2, 1, 0, 0, 0.13 },
+ { "N32", "NH2", "NTE", "N", 2, 3, 3, 0, 0, 0.13 },
+ { "N31", "NH", "NTE", "N", 1, 3, 3, 0, 0, 0.13 },
+ { "N30", "N", "NTE", "N", 0, 3, 3, 0, 0, 0.13 },
+ { "N22", "NH2", "NPI2", "N", 2, 2, 2, 0, 0, 0.12 },
+ { "N21", "NH", "NPI2", "N", 1, 2, 2, 0, 0, 0.12 },
+ { "N20", "N", "NTR2", "N", 0, 2, 2, 0, 0, 0.12 },
+ { "N10", "N", "NDI", "N", 0, 1, 1, 0, 0, 0.12 },
+ { "O31", "OH", "OTE", "O", 1, 2, 3, 0, 0, 0.14 },
+ { "O30", "O", "OTE", "O", 0, 2, 3, 0, 0, 0.14 },
+ { "O20", "O", "OTR4", "O", 0, 1, 2, 0, 0, 0.13 },
+ { "F30", "F", "F", "F", 0, 1, 3, 0, 0, 0.12 },
+ { "Si30", "Si", "Si", "Si", 0, 4, 3, 0, 0, 0.16 },
+ { "P21", "PH", "PTE", "P", 1, 4, 3, 0, 0, 0.17 },
+ { "P20", "P", "PTE", "P", 0, 4, 3, 0, 0, 0.17 },
+ { "S31", "SH", "STE", "S", 1, 4, 3, 0, 0, 0.19 },
+ { "S30", "S", "STE", "S", 0, 4, 3, 0, 0, 0.19 },
+ { "S20", "S", "STR4", "S", 0, 4, 2, 0, 0, 0.18 },
+ { "Cl30", "Cl", "Cl", "Cl", 0, 1, 3, 0, 0, 0.20 },
+ { "Br30", "Br", "Br", "Br", 0, 1, 3, 0, 0, 0.22 },
+ { "I30", "I", "I", "I", 0, 1, 3, 0, 0, 0.24 },
+ { "H", "H", "H", "H", 0, 1, 3, 0, 0, 0.108 },
+ { "He", "He", "He", "He", 0, 0, 0, 0, 0.24, 0 },
+ { "Li+", "Li+", "Li", "Li", 0, 0, 0, 0, 0.032, 0 },
+ { "Be2+", "Be2+", "Be", "Be", 0, 0, 0, 0, 0.0083, 0 },
+ { "B", "B", "B", "B", 0, 0, 0, 0, 0, 0 },
+ { "F-", "F-", "F", "F", 0, 0, 0, 0, 2.467, 0 },
+ { "Ne", "Ne", "Ne", "Ne", 0, 0, 0, 0, 0.44, 0 },
+ { "Na+", "Na+", "Na", "Na", 0, 0, 0, 0, 0.157, 0 },
+ { "Mg2+", "Mg2+", "Mg", "Mg", 0, 0, 0, 0, 0.075, 0 },
+ { "Al3+", "Al3+", "Al", "Al", 0, 0, 0, 0, 0, 0 },
+ { "Cl-", "Cl-", "Cl", "Cl", 0, 0, 0, 0, 5.482, 0 },
+ { "Ar", "Ar", "Ar", "Ar", 0, 0, 0, 0, 1.73, 0 },
+ { "K+", "K+", "K", "K", 0, 0, 0, 0, 0.83, 0 },
+ { "Ca2+", "Ca2+", "Ca", "Ca", 0, 0, 0, 0, 0.49, 0 },
+ { "Br-", "Br-", "Br", "Br", 0, 0, 0, 0, 7.268, 0 },
+ { "Kr", "Kr", "Kr", "Kr", 0, 0, 0, 0, 2.58, 0 },
+ { "Rb+", "Rb+", "Rb", "Rb", 0, 0, 0, 0, 1.37, 0 },
+ { "I-", "I-", "I", "I", 0, 0, 0, 0, 10.275, 0 },
+ { "Xe", "Xe", "Xe", "Xe", 0, 0, 0, 0, 4.15, 0 },
+ { "Cs+", "Cs+", "Cs", "Cs", 0, 0, 0, 0, 2.36, 0 }
+};
+
+t_bosque bosque[eelemNR+1] = {
+ { "H", 1, 0.17, 1.0079 },
+ { "He", 2, 0.0, 4.0026 },
+ { "Li", 3, 0.0, 6.941 },
+ { "Be", 4, 0.0, 9.01218 },
+ { "B", 5, 0.0, 10.81 },
+ { "C", 6, 1.51, 12.011 },
+ { "N", 7, 1.05, 14.00674 },
+ { "O", 8, 0.57, 15.9994 },
+ { "F", 9, 0.22, 18.99984 },
+ { "Ne", 10, 0.0, 20.179 },
+ { "Na", 11, 0.0, 22.98977 },
+ { "Mg", 12, 0.0, 24.305 },
+ { "Al", 13, 0.0, 26.98154 },
+ { "Si", 14, 1.97, 28.086 },
+ { "P", 15, 2.48, 30.97376 },
+ { "S", 16, 2.99, 32.06 },
+ { "Cl", 17, 2.16, 35.453 },
+ { "Ar", 18, 0.0, 39.948 },
+ { "K", 19, 0.0, 39.098 },
+ { "Ca", 20, 0.0, 40.08 },
+ { "Br", 35, 3.29, 79.904 },
+ { "Kr", 36, 0.0, 83.8 },
+ { "Rb", 37, 0.0, 85.478 },
+ { "I", 53, 5.45, 126.9045 },
+ { "Xe", 54, 0.0, 131.3 },
+ { "Cs", 55, 0.0, 132.9054 },
+ { "0", 0, 0.32, 0 }
+};
+
+t_miller miller[emlNR] = {
+ { "H", 1, 0.313, 0.387 },
+ { "F", 9, 1.089, 0.296 },
+ { "Cl", 17, 3.165, 2.315 },
+ { "Br", 35, 5.566, 3.013 },
+ { "I", 53, 8.593, 5.415 },
+ { "CTE", 6, 1.294, 1.061 },
+ { "CTR", 6, 1.433, 1.352 },
+ { "CBR", 6, 1.707, 1.896 },
+ { "CDI", 6, 1.393, 1.283 },
+ { "NTE", 7, 1.373, 0.964 },
+ { "NTR2", 7, 1.262, 1.030 },
+ { "NPI2", 7, 1.220, 1.090 },
+ { "NDI", 7, 1.304, 0.956 },
+ { "OTE", 8, 1.249, 0.637 },
+ { "OTR4", 8, 1.216, 0.569 },
+ { "OPI2", 8, 1.083, 0.274 },
+ { "STE", 16, 3.496, 3.000 },
+ { "STR4", 16, 3.827, 3.729 },
+ { "SPI2", 16, 2.982, 2.700 },
+ { "PTE", 15, 2.485, 1.538 },
+ { "Si", 14, 2.700, 1.97 },
+};
+ #define emlNR asize(miller)
+
+char *lbasis[eqmNR] = {
+ "HF/6-31G", "MP2/aug-cc-pVDZ", "MP2/aug-cc-pVTZ",
+ "B3LYP/aug-cc-pVTZ","B3LYP/d-aug-cc-pVTZ",
+ "Ahc", "Ahp", "Bosque", "TW"
+};
+
+char *ec_name[ecNR] = {
+ "Acid","Alcohol", "Aldehyde", "Alkane", "Alkene", "Alkyne",
+ "Amide", "Amine", "Aromatic",
+ "Ester", "Ether",
+ "Halogen", "Heterocyclic",
+ "Ketone",
+ "Nitril", "Nitro", "Phospho", "Silicate",
+ "Thio", "Other"
+};
+
+int mp_num_prop(t_molprop *mp,char *prop)
+{
+ int i,mnp=0;
+
+ for(i=0; (i<mp->nexperiment); i++) {
+ if (strcasecmp(mp->pname[i],prop) == 0)
+ mnp++;
+ }
+ return mnp;
+}
+
+double mp_get_prop(t_molprop *mp,char *prop,int index)
+{
+ int i,n=0;
+
+ for(i=0; (i<mp->nexperiment); i++) {
+ if (strcasecmp(mp->pname[i],prop) == 0) {
+ if (n == index)
+ return mp->experiment[i];
+ else
+ n++;
+ }
+ }
+ fprintf(stderr,"No %s value for %s (index chosen: %d)\n",
+ prop,mp->molname,index);
+
+ return -1;
+}
+
+char *mp_get_ref_prop(t_molprop *mp,char *prop,int index)
+{
+ int i,n=0;
+
+ for(i=0; (i<mp->nexperiment); i++) {
+ if (strcasecmp(mp->pname[i],prop) == 0) {
+ if (n == index)
+ return mp->reference[i];
+ else
+ n++;
+ }
+ }
+ fprintf(stderr,"No %s value for %s (index chosen: %d)\n",
+ prop,mp->molname,index);
+
+ return NULL;
+}
+
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.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:
- * Good gRace! Old Maple Actually Chews Slate
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdio.h>
-#include <math.h>
-#include "typedefs.h"
-#include "statutil.h"
-#include "copyrite.h"
-#include "gmx_fatal.h"
-#include "xvgr.h"
-#include "pdbio.h"
-#include "macros.h"
-#include "smalloc.h"
-#include "vec.h"
-#include "pbc.h"
-#include "physics.h"
-#include "names.h"
-#include "txtdump.h"
-#include "trnio.h"
-#include "symtab.h"
-#include "confio.h"
-
-real pot(real x,real qq,real c6,real cn,int npow)
-{
- return cn*pow(x,-npow)-c6*pow(x,-6)+qq*ONE_4PI_EPS0/x;
-}
-
-real bhpot(real x,real qq,real A,real B,real C)
-{
- return A*exp(-B*x) - C*pow(x,-6.0);
-}
-
-real dpot(real x,real qq,real c6,real cn,int npow)
-{
- return -(npow*cn*pow(x,-npow-1)-6*c6*pow(x,-7)+qq*ONE_4PI_EPS0/sqr(x));
-}
-
-int main(int argc,char *argv[])
-{
- static char *desc[] = {
- "Plot the potential"
- };
- static real c6=1.0e-3,cn=1.0e-6,qi=0,qj=0,sig=0.3,eps=1,sigfac=0.7;
- static real Abh=1e5,Bbh=32,Cbh=1e-3;
- static int npow=12;
- t_pargs pa[] = {
- { "-c6", FALSE, etREAL, {&c6}, "c6" },
- { "-cn", FALSE, etREAL, {&cn}, "constant for repulsion" },
- { "-pow", FALSE, etINT, {&npow},"power of the repulsion term" },
- { "-sig", FALSE, etREAL, {&sig}, "sig" },
- { "-eps", FALSE, etREAL, {&eps}, "eps" },
- { "-A", FALSE, etREAL, {&Abh}, "Buckingham A" },
- { "-B", FALSE, etREAL, {&Bbh}, "Buckingham B" },
- { "-C", FALSE, etREAL, {&Cbh}, "Buckingham C" },
- { "-qi", FALSE, etREAL, {&qi}, "qi" },
- { "-qj", FALSE, etREAL, {&qj}, "qj" },
- { "-sigfac", FALSE, etREAL, {&sigfac}, "Factor in front of sigma for starting the plot" }
- };
- t_filenm fnm[] = {
- { efXVG, "-o", "potje", ffWRITE }
- };
-#define NFILE asize(fnm)
- static char *legend[] = { "Lennard-Jones", "Buckingham" };
- FILE *fp;
- int i;
- bool bBham;
- real qq,x,oldx,minimum,mval,dp[2],pp[2];
- int cur=0;
-#define next (1-cur)
-
- /* CopyRight(stdout,argv[0]);*/
- parse_common_args(&argc,argv,PCA_CAN_VIEW,
- NFILE,fnm,asize(pa),pa,asize(desc),
- desc,0,NULL);
-
- bBham = (opt2parg_bSet("-A",asize(pa),pa) ||
- opt2parg_bSet("-B",asize(pa),pa) ||
- opt2parg_bSet("-C",asize(pa),pa));
-
- if (bBham) {
- c6 = Cbh;
- sig = pow((6.0/npow)*pow(npow/Bbh,npow-6.0),1.0/(npow-6.0));
- eps = c6/(4*pow(sig,6.0));
- cn = 4*eps*pow(sig,npow);
- }
- else {
- if (opt2parg_bSet("-sig",asize(pa),pa) ||
- opt2parg_bSet("-eps",asize(pa),pa)) {
- c6 = 4*eps*pow(sig,6);
- cn = 4*eps*pow(sig,npow);
- }
- else if (opt2parg_bSet("-c6",asize(pa),pa) ||
- opt2parg_bSet("-cn",asize(pa),pa) ||
- opt2parg_bSet("-pow",asize(pa),pa)) {
- sig = pow(cn/c6,1.0/(npow-6.0));
- eps = 0.25*c6*pow(sig,-6.0);
- }
- else {
- sig = eps = 0;
- }
- printf("c6 = %12.5e, c%d = %12.5e\n",c6,npow,cn);
- printf("sigma = %12.5f, epsilon = %12.5f\n",sig,eps);
-
- minimum = pow(npow/6.0*pow(sig,npow-6.0),1.0/(npow-6));
- printf("Van der Waals minimum at %g, V = %g\n\n",
- minimum,pot(minimum,0,c6,cn,npow));
- printf("Fit of Lennard Jones (%d-6) to Buckingham:\n",npow);
- Bbh = npow/minimum;
- Cbh = c6;
- Abh = 4*eps*pow(sig/minimum,npow)*exp(npow);
- printf("A = %g, B = %g, C = %g\n",Abh,Bbh,Cbh);
- }
- qq = qi*qj;
-
- fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),"Potential","r (nm)","E (kJ/mol)");
- xvgr_legend(fp,asize(legend),legend);
- if (sig == 0)
- sig=0.25;
- minimum = -1;
- mval = 0;
- oldx = 0;
- for(i=0; (i<100); i++) {
- x = sigfac*sig+sig*i*0.02;
- dp[next] = dpot(x,qq,c6,cn,npow);
- fprintf(fp,"%10g %10g %10g\n",x,pot(x,qq,c6,cn,npow),
- bhpot(x,qq,Abh,Bbh,Cbh));
- if (qq != 0) {
- if ((i > 0) && (dp[cur]*dp[next] < 0)) {
- minimum = oldx + dp[cur]*(x-oldx)/(dp[cur]-dp[next]);
- mval = pot(minimum,qq,c6,cn,npow);
- printf("Van der Waals + Coulomb minimum at r = %g (nm). Value = %g (kJ/mol)\n",
- minimum,mval);
- }
- }
- cur = next;
- oldx = x;
-
- }
- fclose(fp);
-
- do_view(ftp2fn(efXVG,NFILE,fnm),NULL);
-
- thanx(stderr);
-
- return 0;
-}
-
-
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
+ * VERSION 3.3.99_development_20071104
* 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,
+ * Copyright (c) 2001-2006, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ * Groningen Machine for Chemical Simulation
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <ctype.h>
+#include <math.h>
#include "maths.h"
#include "macros.h"
#include "copyrite.h"
#include "names.h"
#include "vec.h"
#include "atomprop.h"
-#include "grompp.h"
+#include "xvgr.h"
+#include "viewit.h"
#include "gmx_random.h"
+#include "gstat.h"
#include "x2top_qgen.h"
#include "x2top_eemprops.h"
typedef struct {
char *molname;
- real dip_ref,qtotal,dip_obs;
+ real dip_ref,qtotal,dip_obs,chieq;
t_atoms *atoms;
rvec *x;
matrix box;
} t_mymol;
typedef struct {
- int nmol,nparam;
+ int nmol,nparam,eemtype;
t_mymol *mymol;
- bool bFitJ0,bFitChi0,bFitRadius;
real J0_0,Chi0_0,R_0;
void *eem;
void *atomprop;
mymol->molname = strdup(fn);
mymol->dip_ref = dip;
- mymol->qtotal = 0;
+ mymol->qtotal = 0;
}
-static void print_mol(FILE *fp,t_mymol *mol)
+static void print_mols(FILE *logf,char *xvgfn,int nmol,t_mymol mol[])
{
- int i;
+ FILE *xvgf;
+ double d2=0;
+ real a,b,rms,sigma,error,aver;
+ int i,j,nout;
+ t_lsq lsq;
+
+ xvgf = xvgropen(xvgfn,"Correlation between dipoles",
+ "Experimental","Predicted");
+ fprintf(xvgf,"@ s0 linestyle 0\n");
+ fprintf(xvgf,"@ s0 symbol 1\n");
+ init_lsq(&lsq);
+ for(i=0; (i<nmol); i++) {
+ fprintf(logf,"Molecule %s, Dipole %6.2f, should be %6.2f <chi>: %6.2f\n",
+ mol[i].molname,mol[i].dip_obs,mol[i].dip_ref,mol[i].chieq);
+ fprintf(xvgf,"%10g %10g\n",mol[i].dip_ref,mol[i].dip_obs);
+ add_lsq(&lsq,mol[i].dip_ref,mol[i].dip_obs);
+ d2 += sqr(mol[i].dip_ref-mol[i].dip_obs);
+ fprintf(logf,"Res Atom q\n");
+ for(j=0; (j<mol[i].atoms->nr); j++)
+ fprintf(logf,"%-5s%-5s %8.4f\n",
+ *(mol[i].atoms->resname[mol[i].atoms->atom[j].resnr]),
+ *(mol[i].atoms->atomname[j]),mol[i].atoms->atom[j].q);
+ fprintf(logf,"\n");
+ }
+ fclose(xvgf);
+ get_lsq_ab(&lsq,&a,&b);
+ rms = sqrt(d2/nmol);
+ printf("\nStatistics: fit of %d dipoles Dpred = %.3f Dexp + %3f\n",nmol,a,b);
+ printf("RMSD = %.2f D\n",rms);
+ aver = aver_lsq(&lsq);
+ sigma = rms/aver;
+ nout = 0;
+ fprintf(logf,"Overview of outliers\n");
+ fprintf(logf,"--------------------\n");
+ fprintf(logf,"%-20s %12s %12s %12s\n",
+ "Name","Predicted","Experimental","Deviation");
+ for(i=0; (i<nmol); i++) {
+ if ((mol[i].dip_ref > 0) &&
+ ((mol[i].dip_obs/mol[i].dip_ref-1) > 2*sigma)) {
+ fprintf(logf,"%-20s %12.3f %12.3f %12.3f\n",
+ mol[i].molname,mol[i].dip_obs,mol[i].dip_ref,
+ fabs(mol[i].dip_obs-mol[i].dip_ref));
+ nout ++;
+ }
+ }
+ if (nout)
+ printf("There were %d outliers. See at the very bottom of the log file\n",
+ nout);
+ else
+ fprintf(logf,"None! Well done.\n");
+ do_view(xvgfn,NULL);
+ fclose(logf);
- fprintf(fp,"Molecule %s, Dipole %6.2f, should be %6.2f\n",
- mol->molname,mol->dip_obs,mol->dip_ref);
- fprintf(fp,"Res Atom q\n");
- for(i=0; (i<mol->atoms->nr); i++)
- fprintf(fp,"%-5s%-5s %8.4f\n",
- *(mol->atoms->resname[mol->atoms->atom[i].resnr]),
- *(mol->atoms->atomname[i]),mol->atoms->atom[i].q);
- fprintf(fp,"\n");
+ done_lsq(&lsq);
}
-t_moldip *read_moldip(char *fn,bool bFitJ0,bool bFitChi0,bool bFitRadius,
- real J0_0,real Chi0_0,real R_0)
+t_moldip *read_moldip(FILE *logf,char *fn,char *eem_fn,
+ real J0_0,real Chi0_0,real R_0,int eemtype)
{
char **strings,buf[STRLEN];
int i,nstrings;
gmx_fatal(FARGS,"Error on line %d of %s",i+1,fn);
init_mymol(&(md->mymol[i]),buf,dip);
}
- printf("Read %d sets of molecular coordinates and dipoles\n",nstrings);
- md->eem = read_eemprops(NULL);
- if (md->eem == NULL)
- gmx_fatal(FARGS,"Could not read eemprops.dat");
- md->nparam = eem_getnumprops(md->eem);
+ fprintf(logf,"Read %d sets of molecule coordinates and dipoles\n",nstrings);
+ md->eem = read_eemprops(eem_fn,eemtype);
+ if ((md->eem == NULL) || ((md->nparam = eem_get_numprops(md->eem)) == 0))
+ gmx_fatal(FARGS,"Could not read %s, or file empty",
+ eem_fn ? eem_fn : "eemprops.dat");
+
md->atomprop = get_atomprop();
- md->bFitJ0 = bFitJ0;
- md->bFitChi0 = bFitChi0;
- md->bFitRadius = bFitRadius;
+ md->eemtype = eemtype;
md->J0_0 = J0_0;
md->Chi0_0 = Chi0_0;
md->R_0 = R_0;
-
+ fprintf(logf,"There are %d atom types in the input file %s:\n---\n",
+ md->nparam,
+ eem_fn ? eem_fn : "eemprops.dat");
+ write_eemprops(logf,md->eem);
+ fprintf(logf,"---\n\n");
+
return md;
}
*/
k=0;
for(i=0; (i<md->nparam); i++) {
- if (md->bFitJ0) {
- J0 = gsl_vector_get(v, k++);
- if (J0 <= md->J0_0)
- rms += sqr(J0-md->J0_0);
- }
- else
- J0 = lo_get_j00(md->eem,i,&wj,0);
-
- if (md->bFitChi0) {
- chi0 = gsl_vector_get(v, k++);
- if (chi0 <= md->Chi0_0)
- rms += sqr(chi0-md->Chi0_0);
- }
- else
- chi0 = eem_get_chi0(md->eem,i);
-
- if (md->bFitRadius) {
+ J0 = gsl_vector_get(v, k++);
+ if (J0 <= md->J0_0)
+ rms += sqr(J0-md->J0_0);
+ chi0 = gsl_vector_get(v, k++);
+ if (chi0 <= md->Chi0_0)
+ rms += sqr(chi0-md->Chi0_0);
+ if (md->eemtype != eqgSM1) {
radius = gsl_vector_get(v, k++);
if (radius <= md->R_0)
rms += sqr(radius-md->R_0);
}
for(i=0; (i<md->nmol); i++) {
- generate_charges_sm(debug,md->eem,md->mymol[i].atoms,
- md->mymol[i].x,1e-4,10000,md->atomprop,
- md->mymol[i].qtotal);
+ md->mymol[i].chieq = generate_charges_sm(debug,md->eem,md->mymol[i].atoms,
+ md->mymol[i].x,1e-4,10000,md->atomprop,
+ md->mymol[i].qtotal,md->eemtype);
md->mymol[i].dip_obs = mymol_calc_dip(&(md->mymol[i]));
rms += sqr(md->mymol[i].dip_obs - md->mymol[i].dip_ref);
for(j=0; (j<md->mymol[i].atoms->nr); j++) {
}
static real optimize_moldip(FILE *fp,t_moldip *md,int maxiter,real tol,
- char *outfn,char *logfn)
+ int reinit,real stepsize)
{
- FILE *out;
+ FILE *out,*xvg;
real size,d2,wj;
int iter = 0;
int status = 0;
gsl_multimin_function my_func;
my_func.f = &dipole_function;
- my_func.n = md->nparam*(md->bFitJ0 + md->bFitChi0 + md->bFitRadius);
+ my_func.n = md->nparam*2;
+ if (md->eemtype != eqgSM1)
+ my_func.n += md->nparam;
my_func.params = (void *) md;
/* Starting point */
x = gsl_vector_alloc (my_func.n);
/* Step size, different for each of the parameters */
dx = gsl_vector_alloc (my_func.n);
- k=0;
- for(i=0; (i<md->nparam); i++) {
- if (md->bFitJ0) {
- J00 = lo_get_j00(md->eem,i,&wj,0);
- gsl_vector_set (x, k, J00);
- gsl_vector_set (dx, k++, 0.01*J00);
- }
- if (md->bFitChi0) {
- chi0 = eem_get_chi0(md->eem,i);
- gsl_vector_set (x, k, chi0);
- gsl_vector_set (dx, k++, 0.01*chi0);
- }
- if (md->bFitRadius) {
- radius = eem_get_radius(md->eem,i);
- gsl_vector_set (x, k, radius);
- gsl_vector_set (dx, k++, 0.01*radius);
- }
- }
-
-
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);
-
if (fp)
fprintf(fp,"%5s %12s %12s\n","Iter","Size","RMS");
- do {
+ do {
+ if ((iter == 0) || ((reinit > 0) && (iter % reinit) == 0)) {
+ k=0;
+ for(i=0; (i<md->nparam); i++) {
+ J00 = lo_get_j00(md->eem,i,&wj,0);
+ gsl_vector_set (x, k, J00);
+ gsl_vector_set (dx, k++, stepsize*J00);
+ chi0 = eem_get_chi0(md->eem,i);
+ gsl_vector_set (x, k, chi0);
+ gsl_vector_set (dx, k++, stepsize*chi0);
+
+ if (md->eemtype != eqgSM1) {
+ radius = eem_get_radius(md->eem,i);
+ gsl_vector_set (x, k, radius);
+ gsl_vector_set (dx, k++, stepsize*radius);
+ }
+ }
+ gsl_multimin_fminimizer_set (s, &my_func, x, dx);
+ }
+
iter++;
status = gsl_multimin_fminimizer_iterate (s);
/*if (status == GSL_SUCCESS)
if (fp)
- fprintf(fp,"Minimum found using %s\n",
- gsl_multimin_fminimizer_name(s));
+ fprintf(fp,"Minimum found using %s\n",
+ gsl_multimin_fminimizer_name(s));
*/
if (fp) {
fprintf(fp,"%5d", iter);
}
while (/*(status != GSL_SUCCESS) &&*/ (sqrt(d2) > tol) && (iter < maxiter));
+ gsl_vector_free (x);
+ gsl_vector_free (dx);
gsl_multimin_fminimizer_free (s);
- out = fopen(outfn,"w");
- write_eemprops(out,md->eem);
- fclose(out);
- out = fopen(logfn,"w");
- for(i=0; (i<md->nmol); i++) {
- print_mol(out,&(md->mymol[i]));
- }
- fclose(out);
-
return d2;
}
t_filenm fnm[] = {
{ efDAT, "-f", "moldip", ffREAD },
+ { efDAT, "-d", "eemprops", ffOPTRD },
{ efDAT, "-o", "molprops", ffWRITE },
- { efLOG, "-g", "tune_dip", ffWRITE },
+ { efLOG, "-g", "charges", ffWRITE },
+ { efXVG, "-x", "dipcorr", ffWRITE }
};
#define NFILE asize(fnm)
- static int maxiter=100;
+ static int maxiter=100,reinit=0;
static real tol=1e-3;
- static real J0_0=0,Chi0_0=0,R_0=0;
- static bool bFitJ0=TRUE,bFitChi0=TRUE,bFitRadius=FALSE;
+ static real J0_0=0,Chi0_0=0,R_0=0,step=0.01;
+ static char *qgen[] = { NULL, "SM1", "SM2", "SM3", "SM4", NULL };
t_pargs pa[] = {
{ "-tol", FALSE, etREAL, {&tol},
"Tolerance for convergence in optimization" },
{ "-maxiter",FALSE, etINT, {&maxiter},
"Max number of iterations for optimization" },
- { "-fitj0", FALSE, etBOOL, {&bFitJ0},
- "Optimize dipoles by fitting J00" },
+ { "-reinit", FALSE, etINT, {&reinit},
+ "After this many iterations the search vectors are randomized again. A vlue of 0 means this is never done at all." },
+ { "-qgen", FALSE, etENUM, {qgen},
+ "Algorithm used for charge generation" },
{ "-j0", FALSE, etREAL, {&J0_0},
"Minimum value that J0 can obtain in fitting" },
- { "-fitChi0", FALSE, etBOOL, {&bFitChi0},
- "Optimize dipoles by fiting Chi0" },
{ "-chi0", FALSE, etREAL, {&Chi0_0},
"Minimum value that Chi0 can obtain in fitting" },
- { "-fitradius", FALSE, etBOOL, {&bFitRadius},
- "Optimize dipole by fitting Radius" },
{ "-r0", FALSE, etREAL, {&R_0},
- "Minimum value that Radius can obtain in fitting" }
+ "Minimum value that Radius can obtain in fitting" },
+ { "-step", FALSE, etREAL, {&step},
+ "Step size in parameter optimization. Is used as a fraction of the starting value, should be less than 10%. At each reinit step the step size is updated." }
};
t_moldip *md;
-
+ int eemtype;
+ FILE *logf,*out;
+
CopyRight(stdout,argv[0]);
- parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
+ parse_common_args(&argc,argv,PCA_CAN_VIEW,NFILE,fnm,asize(pa),pa,
asize(desc),desc,0,NULL);
- if (!bFitJ0 && !bFitChi0 && !bFitRadius)
- gmx_fatal(FARGS,"Nothing to fit to!");
-
- md = read_moldip(opt2fn("-f",NFILE,fnm),bFitJ0,bFitChi0,bFitRadius,
- J0_0,Chi0_0,R_0);
-
- (void) optimize_moldip(stdout,md,maxiter,tol,
- opt2fn("-o",NFILE,fnm),
- opt2fn("-g",NFILE,fnm));
+ eemtype = eqgSM1;
+ if (qgen[0]) {
+ eemtype = name2eemtype(qgen[0]);
+ if (eemtype == -1)
+ eemtype = eqgSM1;
+ }
+ if (eemtype > eqgSM3)
+ gmx_fatal(FARGS,"Only models SM1 and SM2 implemented so far");
+ logf = fopen(opt2fn("-g",NFILE,fnm),"w");
+ md = read_moldip(logf,opt2fn("-f",NFILE,fnm),opt2fn_null("-d",NFILE,fnm),
+ J0_0,Chi0_0,R_0,eemtype);
+
+ (void) optimize_moldip(stdout,md,maxiter,tol,reinit,step);
+
+ print_mols(logf,opt2fn("-x",NFILE,fnm),md->nmol,md->mymol);
+
+ fclose(logf);
+
+ out = fopen(opt2fn("-o",NFILE,fnm),"w");
+ write_eemprops(out,md->eem);
+ fclose(out);
+
thanx(stderr);
return 0;
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.3.99_development_20071104
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2006, 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:
+ * Groningen Machine for Chemical Simulation
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "maths.h"
+#include "tune_pol.h"
+#include "vec.h"
+#include "x2top_matrix.h"
+
+static int tabnum = 0;
+static int toler = 5; /* Tolerance in % */
+static int tolerb = 50; /* Tolerance in % */
+static int my_order[eelemNR] = {
+ eelemC,
+ eelemH, eelemHe,
+ eelemLi, eelemBe, eelemB, eelemN, eelemO, eelemF, eelemNe,
+ eelemNa, eelemMg, eelemAl, eelemSi, eelemP, eelemS, eelemCl, eelemAr,
+ eelemK, eelemCa, eelemBr, eelemKr, eelemRb,
+ eelemI, eelemXe, eelemCs
+};
+
+/*static char *sbasis[eqmNR] = {
+ "6-31G","pVDZ","pVTZ","pVTZ","d-pVTZ", "Ahc", "Ahp", "BS", "FP"
+ };*/
+
+static char *sbasis[eqmNR] = {
+ "6-31G","aug-DZ","aug-TZ","aug-TZ","d-aug-TZ", "Ahc", "Ahp", "BS", "FP"
+};
+
+static void check_formula(int np,t_molprop pd[])
+{
+ int i,j,jj;
+ char formula[1280],number[32];
+
+ for(i=0; (i<np); i++) {
+ formula[0] = '\0';
+ for(j=0; (j<eelemNR); j++) {
+ jj = my_order[j];
+ if (pd[i].elem_comp[jj] > 0) {
+ strcat(formula,bosque[jj].name);
+ if (pd[i].elem_comp[jj] > 1) {
+ sprintf(number,"%d",pd[i].elem_comp[jj]);
+ strcat(formula,number);
+ }
+ }
+ }
+ if (0 && (strcasecmp(pd[i].formula,formula) != 0) )
+ fprintf(stderr,"Formula '%s' does match '%s' based on composition for %s.\n",
+ pd[i].formula,formula,pd[i].molname);
+ }
+}
+
+static void table_number(FILE *fp)
+{
+ fprintf(fp,"\\setcounter{table}{%d}\n",tabnum);
+}
+
+static void table1_header(FILE *fp,int caption)
+{
+ table_number(fp);
+ fprintf(fp,"\\begin{table}[H]\n\\centering\n");
+ if (caption)
+ fprintf(fp,"\\caption{Definition of fragment-based atom types. The naming is such that the first digit represents the hybridization state and the second digit the number of hydrogen atoms attached. The number of occurences of the fragments N is given along with the fragment polarizabilities (FP). The columns Ahc and Ahp contain group polarizabilites computed using Miller's equation~\\protect\\cite{Miller1979a} and parameters~\\protect\\cite{Miller90a} and Kang and Jhon's method~\\cite{Kang1982a} with the parametrization of Miller~\\protect\\cite{Miller90a} respectively. The column BS contains the equivalent using the polarizabilities of Bosque and Sales~\\protect\\cite{Bosque2002a}.}\n\\label{fragments}\n");
+ else
+ fprintf(fp,"\\caption{Continued}");
+ fprintf(fp,"\\begin{tabular}{cccccccc}\n\\hline\n");
+ fprintf(fp,"Group & Name & Hybridiz. & \\multicolumn{1}{c}{N} & \\multicolumn{4}{c}{Polarizability} \\\\\n");
+ fprintf(fp,"& & & & Ahc & Ahp & BS & FP \\\\\n");
+ fprintf(fp,"\\hline\n");
+}
+
+static void table1_end(FILE *fp)
+{
+ fprintf(fp,"\\hline\n");
+ fprintf(fp,"\\end{tabular}\n\\end{table}\n\n");
+}
+
+static int miller_index(char *mil)
+{
+ int imil;
+
+ for(imil=0; (imil<emlNR); imil++)
+ if (strcmp(mil,miller[imil].name) == 0)
+ break;
+ if (imil == emlNR) {
+ fprintf(stderr,"No such Miller atom %s\n",mil);
+ exit(1);
+ }
+ return imil;
+}
+
+static int bosque_index(char *bos)
+{
+ int ielem;
+
+ for(ielem=0; (ielem<eelemNR); ielem++)
+ if (strcmp(bos,bosque[ielem].name) == 0)
+ break;
+ if (ielem == eelemNR) {
+ fprintf(stderr,"No such element %s",bos);
+ exit(1);
+ }
+ return ielem;
+}
+
+static void dump_table1(FILE *fp,int np,t_molprop pd[])
+{
+ int i,j,nfit,nh,ielem,imil;
+ char hybrid[8];
+ double ahc,bos,ahp;
+
+ table1_header(fp,1);
+ strcpy(hybrid,"SP ");
+
+ for(i=0; (i<eatNR); i++) {
+ int len = strlen(spoel[i].name);
+ if (len < 3) {
+ fprintf(stderr,"No such fragment %s",spoel[i].name);
+ }
+ hybrid[2] = spoel[i].name[len-2];
+ nh = spoel[i].nh;
+ ielem = bosque_index(spoel[i].bosque);
+ imil = miller_index(spoel[i].miller);
+ bos = bosque[ielem].value + nh*bosque[eelemH].value;
+ nfit = 0;
+ for(j=0; (j<np); j++) {
+ if ((mp_num_polar(&(pd[j])) > 0) && (mp_get_polar(&(pd[j])) > 0))
+ nfit += pd[j].frag_comp[i];
+ }
+ ahc = (4.0/(nh+miller[imil].nelec))*sqr(miller[imil].tau_ahc +
+ nh*miller[emlH].tau_ahc);
+ ahp = (nh*miller[emlH].alpha_ahp + miller[imil].alpha_ahp);
+ fprintf(fp,"%s & %s & %s & %d & %8.2f & %8.2f & %8.2f & %8.2f \\\\\n",
+ spoel[i].atom,spoel[i].name,hybrid,nfit,
+ ahc,ahp,bos,spoel[i].all);
+ if ((((i+1) % 16) == 0) && (i != eatNR-1)) {
+ table1_end(fp);
+ table1_header(fp,0);
+ }
+ }
+ table1_end(fp);
+}
+
+static void table2_header(FILE *fp,int caption,int bTrain,int bQM)
+{
+ int i;
+ char *ptr;
+
+ ptr = bTrain ? "training" : "test";
+ table_number(fp);
+ fprintf(fp,"\\begin{%stable}[H]\n",bQM ? "sideways" : "");
+ if (caption)
+ fprintf(fp,"\\caption{Molecules with experimental polarizabilities as well as %sfragment or atom based polarizabilities of Miller {\\em et al.}~\\protect\\cite{Miller1979a,Miller90a} (Ahc, Ahp), Bosque and Sales~\\protect\\cite{Bosque2002a} (BS) or fragment polarizabilities (FP, this work). %sCalculated numbers that are more than %d\\%% off are printed in bold.%s A (*) behind the experimental value indicates that these were derived in this work. %s}\n\\label{%s}",
+ bQM ? "quantum chemical and " : "",
+ bQM ? "Basissets used are aug-cc-pVXZ and d-aug-cc-pVXZ. " : "",
+ toler,
+ bQM ? " Numbering is identical to that in Supplementary material, Table 1." : "",
+ bQM ? "$\\dagger$: for this calculation no convergence was achieved and hence no polarizability could be calculated." : "",
+ bQM ? "qm" : "empirical");
+ else
+ fprintf(fp,"\\caption{Continued}\n");
+ if (bQM) {
+ fprintf(fp,"\\begin{tabular}{lcccccccccc}\n\\hline\n");
+ fprintf(fp,"Molecule & Experiment & HF & \\multicolumn{2}{c}{MP2} & \\multicolumn{2}{c}{B3LYP} & \\multicolumn{4}{c}{Empirical}\\\\\n");
+ fprintf(fp," & ");
+ for(i=0; (i<eqmNR); i++)
+ fprintf(fp,"& %s ",sbasis[i]);
+ fprintf(fp,"\\\\\n\\hline\n");
+ }
+ else {
+ fprintf(fp,"\\begin{tabular}{lccccc}\n\\hline\n");
+ fprintf(fp,"Molecule & Experiment ");
+ for(i=NQUANT; (i<eqmNR); i++)
+ fprintf(fp,"& %s ",sbasis[i]);
+ fprintf(fp,"\\\\\n\\hline\n");
+ }
+}
+
+static void table2_end(FILE *fp,int bQM)
+{
+ fprintf(fp,"\\hline\n\\end{tabular}\n\\end{%stable}\n",
+ bQM ? "sideways" : "");
+}
+
+static void dump_table2(FILE *fp,int bVerbose,int bTrain,int bQM,
+ int np,t_molprop pd[],int bDS)
+{
+ int i,j,j0,k,header,iline,iqm,caption=1,maxline;
+ char buf[32],buf2[32],*ref;
+ double val;
+
+ table2_header(fp,caption,bTrain,bQM);
+ header = 1;
+ iline = 0;
+ j0 = bQM ? 0 : NQUANT;
+ for(i=0; (i<np); i++) {
+ iqm = 0;
+ for(k=0; (k<NQUANT); k++)
+ if (pd[i].qm[k] > 0)
+ iqm++;
+ if (((bTrain && (pd[i].weight > 0)) ||
+ (!bTrain && (pd[i].weight == 0))) &&
+ (mp_num_polar(&(pd[i])) > 0) &&
+ ((bQM && (iqm > 0)) || !bQM)) {
+ fprintf(fp,"%3d. %-15s",pd[i].nr,pd[i].molname);
+ header = 0;
+ fprintf(fp," &");
+ for(k=0; (k<mp_num_polar(&(pd[i]))); k++,iline++) {
+ if (k > 0)
+ fprintf(fp," &");
+ val = mp_get_prop(&(pd[i]),"Polarizability",k);
+ ref = mp_get_ref_prop(&(pd[i]),"Polarizability",k);
+ if (val > 0) {
+ fprintf(fp," %8.3f",val);
+ if (strcmp(ref,"Spoel2007a") == 0)
+ fprintf(fp," (*)");
+ else
+ fprintf(fp,"~\\cite{%s} ",ref ? ref : "XXX");
+ }
+ if (k == 0) {
+ for(j=j0; (j<eqmNR); j++)
+ if (pd[i].qm[j] > 0) {
+ if ((fabs(pd[i].qm[j] - val) > (val*toler*0.01))) {
+ if (fabs(pd[i].qm[j] - val) > (val*tolerb*0.01))
+ fprintf(fp,"& {\\LARGE\\bf %8.2f} ",pd[i].qm[j]);
+ else
+ fprintf(fp,"& {\\bf %8.2f} ",pd[i].qm[j]);
+ }
+ else
+ fprintf(fp,"& %8.2f ",pd[i].qm[j]);
+ }
+ else
+ fprintf(fp,"& $\\dagger$");
+ fprintf(fp,"\\\\\n");
+ }
+ else {
+ if (bQM)
+ fprintf(fp,"&&&&&&&&&\\\\\n");
+ else
+ fprintf(fp,"&&&&\\\\\n");
+ }
+ }
+ maxline = 28;
+ if (!bQM)
+ maxline += 8;
+ if (caption == 1)
+ maxline -= 5;
+ if (bDS)
+ maxline = maxline*0.6;
+
+ if ((iline >= maxline) && (i+1<np)) {
+ caption = 0;
+ table2_end(fp,bQM);
+ table2_header(fp,caption,bTrain,bQM);
+ header = 1;
+ iline = 0;
+ }
+ }
+ }
+ table2_end(fp,bQM);
+}
+
+static void dump_table3(FILE *fp,int np,t_molprop pd[],int ntot)
+{
+ int i,j,k,ns,n,nn[eqmNR];
+ double aa[eqmNR],bb[eqmNR],R[eqmNR],rms,diff,ssx,ssy,ssxy,ssxx,ssyy;
+ double ra[eqmNR],pp[eqmNR],ratio[eqmNR];
+ double val;
+
+ table_number(fp);
+ fprintf(fp,"\\begin{table}[H]\n\\centering\n");
+ fprintf(fp,"\\caption{Performance of the different methods for predicting molecular polarizability for different subsets of the test molecules. RMSD from experimental values (\\AA$^3$), in brackets the number of molecules in this particular subset. For the empirical methods (Ahc,Ahp~\\cite{Miller90a}, BS~\\cite{Bosque2002a} and FP (this work) all molecules were taken into account. At the bottom the \\%% difference, the correlation coefficient R, the regression coefficient a and the intercept b are given.}\n\\label{frag_rmsd}\n");
+ fprintf(fp,"\\begin{tabular}{lcccccccccc}\n\\hline\n");
+ fprintf(fp,"Method & HF & \\multicolumn{2}{c}{MP2} & \\multicolumn{2}{c}{B3LYP} & \\multicolumn{4}{c}{Empirical}\\\\\n");
+
+ for(i=0; (i<eqmNR); i++)
+ fprintf(fp,"& %s ",sbasis[i]);
+ fprintf(fp,"\\\\\n\\hline\n");
+
+ for(i=0; (i<ecNR); i++) {
+ n = 0;
+ for(j=0; (j<np); j++)
+ if (pd[j].category[i] > 0)
+ n++;
+ fprintf(fp," %s (%d) ",ec_name[i],n);
+ for(k=0; (k<eqmNR); k++) {
+ rms = 0;
+ n = 0;
+ for(j=0; (j<np); j++) {
+ if (mp_num_polar(&(pd[j])) > 0) {
+ val = mp_get_polar(&(pd[j]));
+ if ((val > 0) && (pd[j].qm[k] > 0) && (pd[j].category[i] > 0)) {
+ rms += sqr(val - pd[j].qm[k]);
+ n++;
+ }
+ }
+ }
+ if (n > 0) {
+ if (k >= NQUANT)
+ fprintf(fp,"& %8.2f",sqrt(rms/n));
+ else
+ fprintf(fp,"& %8.2f(%d)",sqrt(rms/n),n);
+ }
+ else
+ fprintf(fp,"& -");
+ }
+ fprintf(fp,"\\\\\n");
+ }
+ for(k=0; (k<eqmNR); k++) {
+ ra[k] = 0;
+ pp[k] = 0;
+ nn[k] = 0;
+ ratio[k] = 0;
+ for(j=0; (j<np); j++) {
+ if (mp_num_polar(&(pd[j])) > 0) {
+ val = mp_get_polar(&(pd[j]));
+ if ((val > 0) && (pd[j].qm[k] > 0)) {
+ diff = val - pd[j].qm[k];
+ ra[k] += sqr(diff);
+ pp[k] += 100*diff/val;
+ ratio[k] += pd[j].qm[k]/val;
+ nn[k]++;
+ }
+ }
+ }
+ }
+ fprintf(fp,"All(%d)",nn[eqmNR-1]);
+ for(k=0; (k<eqmNR); k++) {
+ if (k < NQUANT)
+ fprintf(fp,"& %8.2f(%d)",sqrt(ra[k]/nn[k]),nn[k]);
+ else
+ fprintf(fp,"& %8.2f",sqrt(ra[k]/nn[k]));
+ }
+ fprintf(fp,"\\\\\n\\hline\n");
+ fprintf(fp,"\\%% Diff ");
+ for(k=0; (k<eqmNR); k++) {
+ if (nn[k] > 0) {
+ fprintf(fp,"& %8.2f",pp[k]/nn[k]);
+ }
+ else
+ fprintf(fp,"& -");
+ }
+ fprintf(fp,"\\\\\n");
+ for(k=0; (k<eqmNR); k++) {
+ ssx = ssy = ssxx = ssxy = ssyy = 0;
+ ns = 0;
+ for(j=0; (j<np); j++) {
+ if (mp_num_polar(&(pd[j])) > 0) {
+ val = mp_get_polar(&(pd[j]));
+ if ((val > 0) && (pd[j].qm[k] > 0)) {
+ ssx += val;
+ ssy += pd[j].qm[k];
+ ssxx += sqr(val);
+ ssyy += sqr(pd[j].qm[k]);
+ ssxy += val*pd[j].qm[k];
+ ns++;
+ }
+ }
+ }
+ ssx /= ns;
+ ssy /= ns;
+ ssxy -= ns*ssx*ssy;
+ ssxx -= ns*sqr(ssx);
+ ssyy -= ns*sqr(ssy);
+ R[k] = 100*(ssxy*ssxy)/(ssxx*ssyy);
+ bb[k] = ssxy/ssxx;
+ aa[k] = ssy-bb[k]*ssx;
+ }
+ fprintf(fp,"R (\\%%) ");
+ for(k=0; (k<eqmNR); k++) {
+ fprintf(fp,"& %8.1f",R[k]);
+ }
+ fprintf(fp,"\\\\\n");
+ fprintf(fp,"a ");
+ for(k=0; (k<eqmNR); k++) {
+ fprintf(fp,"& %8.2f",bb[k]);
+ }
+ fprintf(fp,"\\\\\n");
+ fprintf(fp,"b ");
+ for(k=0; (k<eqmNR); k++) {
+ fprintf(fp,"& %8.2f",aa[k]);
+ }
+ fprintf(fp,"\\\\\n");
+ fprintf(fp,"\\hline\n\\end{tabular}\n\\end{table}\n");
+}
+
+static void table4_header(FILE *fp,int caption)
+{
+ table_number(fp);
+ fprintf(fp,"\\begin{sidewaystable}[H]\n\\centering\n");
+ if (caption)
+ fprintf(fp,"\\caption{Decomposition of molecules into fragments or Miller atoms~\\protect\\cite{Miller1979a,Miller90a}}\n\\label{frag_defs}\n");
+ else
+ fprintf(fp,"\\caption{Continued}\n");
+ fprintf(fp,"\\begin{tabular}{llll}\n\\hline\n");
+ fprintf(fp,"Molecule & Formula & Fragments & Miller \\\\\n");
+ fprintf(fp,"\\hline\n");
+}
+
+static void dump_table4(FILE *fp,int np,t_molprop pd[])
+{
+ int i,j,k,iline;
+ char buf[32],buf2[32];
+
+ table4_header(fp,1);
+ iline = 0;
+ for(i=0; (i<np); i++) {
+ if ((i == 0) ||
+ ((i > 0) && strcmp(pd[i].molname,pd[i-1].molname) != 0)) {
+ fprintf(fp,"%3d. %s & %s &",++iline,pd[i].molname,pd[i].formula);
+ for(j=0; (j<eatNR); j++)
+ if (pd[i].frag_comp[j] > 0)
+ fprintf(fp,"%s$_{%d}$ ",spoel[j].name,pd[i].frag_comp[j]);
+ fprintf(fp," &");
+ for(j=0; (j<emlNR); j++)
+ if (pd[i].emil_comp[j] > 0)
+ fprintf(fp,"%s$_{%d}$ ",miller[j].name,pd[i].emil_comp[j]);
+ fprintf(fp," \\\\\n");
+ if (((iline % 28) == 0) && (i+1<np)) {
+ fprintf(fp,"\\hline\n\\end{tabular}\n\\end{sidewaystable}\n");
+ table4_header(fp,0);
+ }
+ }
+ }
+ fprintf(fp,"\\hline\n\\end{tabular}\n\\end{sidewaystable}\n");
+}
+
+static void calc_frag_miller(int bTrain,int np,t_molprop pd[])
+{
+ int j,k,Nelec;
+ double ahc,ahp,bos,spol;
+
+ for(j=0; (j<np); j++) {
+ spol = 0;
+ for(k=0; (k<eatNR); k++)
+ spol += pd[j].frag_comp[k] * (bTrain ? spoel[k].train : spoel[k].all);
+ pd[j].qm[eqmSpoel] = spol;
+ bos = bosque[eelemNR].value;
+ for(k=0; (k<eelemNR); k++)
+ bos += bosque[k].value*pd[j].elem_comp[k];
+ if ((pd[j].qm[eqmBosque] != 0) &&
+ (fabs(pd[j].qm[eqmBosque] - bos) >= 0.01*toler*bos)) {
+ fprintf(stderr,"Possible error in Bosque composition for %s. Computed %g iso %g\n",pd[j].molname,bos,pd[j].qm[eqmBosque]);
+ }
+ pd[j].qm[eqmBosque] = bos;
+ ahc = 0;
+ Nelec = 0;
+ ahp = 0;
+ for(k=0; (k<emlNR); k++) {
+ ahc += miller[k].tau_ahc*pd[j].emil_comp[k];
+ ahp += miller[k].alpha_ahp*pd[j].emil_comp[k];
+ Nelec += miller[k].nelec*pd[j].emil_comp[k];
+ }
+ pd[j].qm[eqmMiller] = 4*sqr(ahc)/Nelec;
+ pd[j].qm[eqmKang] = ahp;
+ }
+}
+
+static void calc_rmsd(int bTrain,double fit[],double test[],int np,
+ t_molprop pd[])
+{
+ int i,j,k,nfit[eqmNR],ntest[eqmNR],Nelec;
+ double eval,qval;
+ double fpol,apol;
+
+ calc_frag_miller(bTrain,np,pd);
+
+ for(i=0; (i<eqmNR); i++) {
+ nfit[i] = 0;
+ ntest[i] = 0;
+ fit[i] = 0;
+ test[i] = 0;
+ for(j=0; (j<np); j++) {
+ if (mp_num_polar(&(pd[j])) > 0) {
+ eval = mp_get_polar(&(pd[j]));
+ if (eval > 0) {
+ qval = pd[j].qm[i];
+
+ if (qval > 0) {
+ if (pd[j].weight > 0) {
+ nfit[i]++;
+ fit[i] += sqr(qval-eval);
+ }
+ else {
+ ntest[i]++;
+ test[i] += sqr(qval-eval);
+ }
+ }
+ }
+ }
+ }
+ if (nfit[i] > 0)
+ fit[i] = sqrt(fit[i]/nfit[i]);
+ if (ntest[i] > 0)
+ test[i] = sqrt(test[i]/ntest[i]);
+ }
+}
+
+static void decompose_frag(FILE *fp,int bTrain,int np,t_molprop pd[])
+{
+ double *x,*atx;
+ double **a,**at,**ata,fpp;
+ int i,j,n,m=eatNR,nn;
+ int test[eatNR];
+
+ for(i=0; (i<eatNR); i++) {
+ test[i] = 0;
+ for(j=0; (j<np); j++)
+ if (pd[j].weight > 0)
+ test[i] += pd[j].frag_comp[i];
+ if (test[i] == 0) {
+ fprintf(stderr,"You have no molecules with group %s\n",
+ spoel[i].name);
+ exit(1);
+ }
+ }
+ x = calloc(np,sizeof(x[0]));
+ for(i=n=0; (i<np); i++) {
+ if ((pd[i].weight > 0) &&
+ (mp_num_polar(&(pd[i])) > 0)) {
+ x[n] = mp_get_polar(&(pd[i]));
+ n++;
+ }
+ }
+ a = alloc_matrix(n,m);
+ at = alloc_matrix(m,n);
+ ata = alloc_matrix(m,m);
+ for(i=nn=0; (i<np); i++) {
+ if ((pd[i].weight > 0) &&
+ (mp_num_polar(&(pd[i])) > 0)){
+ for(j=0; (j<m); j++) {
+ a[nn][j] = pd[i].frag_comp[j];
+ at[j][nn] = pd[i].frag_comp[j];
+ }
+ nn++;
+ }
+ }
+ if (n != nn) {
+ fprintf(stderr,"Consistency error %s, line %d\n",__FILE__,__LINE__);
+ exit(1);
+ }
+ matrix_multiply(fp,n,m,a,at,ata);
+ matrix_invert(fp,m,ata);
+ atx = calloc(m,sizeof(atx[0]));
+ for(i=0; (i<m); i++)
+ for(j=0; (j<n); j++)
+ atx[i] += at[i][j]*x[j];
+
+ for(i=0; (i<m); i++) {
+ fpp = 0;
+ for(j=0; (j<m); j++)
+ fpp += ata[i][j]*atx[j];
+ if (bTrain)
+ spoel[i].train = fpp;
+ else
+ spoel[i].all = fpp;
+ }
+}
+
+static void write_xvg(int np,t_molprop pd[])
+{
+ char *fn[2] = { "qm.xvg", "empir.xvg" };
+ FILE *fp;
+ int i,j,n,i0,i1;
+ double val;
+
+ for(n=0; (n<2); n++) {
+ fp = fopen(fn[n],"w");
+ if (n == 0) {
+ i0 = 0;
+ i1 = NQUANT;
+ val = 19.5;
+ }
+ else {
+ i0 = NQUANT;
+ i1 = eqmNR;
+ val = 51.5;
+ }
+ for(i=i0; (i<i1); i++) {
+ fprintf(fp,"@type xy\n");
+ for(j=0; (j<np); j++) {
+ if ((mp_num_polar(&(pd[j])) > 0) &&
+ (mp_get_polar(&(pd[j])) > 0) &&
+ (pd[j].qm[i] > 0))
+ fprintf(fp,"%8.3f %8.3f\n",mp_get_polar(&(pd[j])),pd[j].qm[i]);
+ }
+ fprintf(fp,"&\n");
+ }
+ fprintf(fp,"@type xy\n");
+ fprintf(fp,"1.5 1.5\n%6.1f %6.1f\n",val,val);
+ fprintf(fp,"&\n");
+ fclose(fp);
+ }
+}
+
+static void write_csv(int np,t_molprop pd[])
+{
+ FILE *fp,*gp;
+ int i,j;
+
+ fp = fopen("matrix.csv","w");
+ gp = fopen("vector.csv","w");
+ for(i=0; (i<np); i++) {
+ if ((pd[i].weight > 0) &&
+ (mp_num_polar(&(pd[i])) > 0) &&
+ (mp_get_polar(&(pd[i])) > 0)) {
+ fprintf(gp,"%10.3f\n",mp_get_polar(&(pd[i])));
+ for(j=0; (j<eatNR-1); j++)
+ fprintf(fp," %5d,",pd[i].frag_comp[j]);
+ fprintf(fp," %5d\n",pd[i].frag_comp[eatNR-1]);
+ }
+ }
+ fclose(fp);
+ fclose(gp);
+}
+
+static int comp_pd_name(const void *a,const void *b)
+{
+ t_molprop *pa,*pb;
+ pa = (t_molprop *)a;
+ pb = (t_molprop *)b;
+
+ return strcasecmp(pa->molname,pb->molname);
+}
+
+static int comp_pd_elem(const void *a,const void *b)
+{
+ t_molprop *pa,*pb;
+ pa = (t_molprop *)a;
+ pb = (t_molprop *)b;
+
+ return strcmp(pa->formula,pb->formula);
+}
+
+static int comp_pd_elem2(const void *a,const void *b)
+{
+ int i,oo=0,nn=0;
+ t_molprop *pa,*pb;
+ pa = (t_molprop *)a;
+ pb = (t_molprop *)b;
+
+ for(i=2; ((i<eelemNR) && (nn == 0)); i++)
+ if ((pa->elem_comp[i] == 0) && (pb->elem_comp[i] > 0)) {
+ nn = -1;
+ break;
+ }
+ else if ((pb->elem_comp[i] == 0) && (pa->elem_comp[i] > 0)) {
+ nn = 1;
+ break;
+ }
+
+ if (nn != 0)
+ return nn;
+
+ for(i=0; (i<eelemNR); i++)
+ if ((oo = (pa->elem_comp[my_order[i]] - pb->elem_comp[my_order[i]])) != 0)
+ return oo;
+
+ return strcasecmp(pa->molname,pb->molname);
+}
+
+static void dump_n2t(char *fn)
+{
+ FILE *fp;
+ int i,j,ielem;
+ double q = 0, m = 0;
+
+ fp = fopen(fn,"w");
+ for(i=0; (i<eatNR+eatExtra); i++) {
+ ielem = bosque_index(spoel[i].bosque);
+ fprintf(fp,"%3s %4s %9g %9g %1d",
+ spoel[i].bosque,spoel[i].name,
+ spoel[i].all,bosque[ielem].mass,
+ spoel[i].maxbond);
+ for(j=0; (j<spoel[i].nh); j++)
+ fprintf(fp," H %.3f",spoel[eatNR].blen);
+ for( ; (j<spoel[i].maxbond); j++)
+ fprintf(fp," * %.3f",spoel[i].blen);
+ fprintf(fp,"\n");
+ }
+
+ fclose(fp);
+}
+
+static int read_dipoles(char *fn)
+{
+ FILE *fp;
+
+ fp = fopen(fn,"r");
+
+ fclose(fp);
+
+ return 0;
+}
+
+int main(int argc,char *argv[])
+{
+ FILE *fp;
+ int i,j,np,nspoel,nbosque,ntot,nqm,nqqm,nqmtot=0;
+ double fit1[eqmNR],test1[eqmNR];
+ double fit2[eqmNR],test2[eqmNR];
+ double *w;
+ t_molprop *mp=NULL;
+
+ if (argc < 2)
+ fatal_error("Please give names of database file!","");
+ mp = merge_xml(argc,argv,NULL,NULL,&np);
+
+ check_formula(np,mp);
+ snew(w,np);
+ qsort(mp,np,sizeof(mp[0]),comp_pd_elem2);
+
+ for(i=0; (i<np); i++) {
+ mp[i].nr = 1+i;
+ w[i] = mp[i].weight;
+ if ((mp[i].weight > 0) &&
+ (mp_num_polar(&(mp[i])) > 0) &&
+ (mp_get_polar(&(mp[i])) == 0)) {
+ fprintf(stderr,"Inconsistency for %s: weight = %g, experiment = 0\n",
+ mp[i].molname,mp[i].weight);
+ exit(1);
+ }
+ }
+
+ nbosque = nspoel = nqm = ntot = 0;
+ for(i=0; (i<np); i++) {
+ if ((mp_num_polar(&(mp[i])) > 0) && (mp_get_polar(&(mp[i])) > 0)) {
+ char *ref = mp_get_ref_polar(&(mp[i]));
+ mp[i].weight = 1;
+ if (strcasecmp(ref,"Spoel2007a") == 0)
+ nspoel++;
+ else if (strcasecmp(ref,"Bosque2002a") == 0)
+ nbosque++;
+ nqqm = 0;
+ for(j=0; (j<NQUANT); j++)
+ if (mp[i].qm[j] > 0)
+ nqqm++;
+ if (nqqm)
+ nqm++;
+ nqmtot += nqqm;
+ ntot++;
+ }
+ else {
+ /*fprintf(stderr,"No experimental data for %s\n",mp[i].molname);*/
+ }
+ if ((i > 0) && (strcasecmp(mp[i].molname,mp[i-1].molname) == 0))
+ fprintf(stderr,"Double entry %s\n",mp[i].molname);
+ }
+ printf("--------------------------------------------------\n");
+ printf(" Some Statistics\n");
+ printf("There are %d entries with experimental data\n",ntot);
+ printf("There are %d experiments with Spoel2007a as reference\n",nspoel);
+ printf("There are %d experiments with Bosque2002a as reference\n",nbosque);
+ printf("There are %d entries with quantum chemistry results, %d missing\n",
+ nqm,nqm*5-nqmtot);
+ printf("There are %d Spoel atom types\n",eatNR);
+ printf("--------------------------------------------------\n");
+
+ decompose_frag(NULL,0,np,mp);
+ calc_rmsd(0,fit2,test2,np,mp);
+
+ if ((fp=fopen("molecules.tex","w")) == NULL)
+ exit(1);
+
+ dump_table1(fp,np,mp);
+ tabnum++;
+ dump_table2(fp,0,1,1,np,mp,1);
+ tabnum++;
+ dump_table3(fp,np,mp,ntot);
+ tabnum = 0;
+ fclose(fp);
+
+ if ((fp=fopen("molecules-sup.tex","w")) == NULL)
+ exit(1);
+ dump_table2(fp,0,1,0,np,mp,0);
+ tabnum++;
+ dump_table4(fp,np,mp);
+ fclose(fp);
+ dump_n2t("ffsm.n2t");
+ write_xvg(np,mp);
+
+ return 0;
+}
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.3.99_development_20071104
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2006, 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:
+ * Groningen Machine for Chemical Simulation
+ */
+
+#ifndef _tune_pol_h
+#define _tune_pol_h
+
+enum {
+ eelemH, eelemHe,
+ eelemLi, eelemBe, eelemB, eelemC, eelemN, eelemO, eelemF, eelemNe,
+ eelemNa, eelemMg, eelemAl, eelemSi, eelemP, eelemS, eelemCl, eelemAr,
+ eelemK, eelemCa, eelemBr, eelemKr,
+ eelemRb, eelemI, eelemXe, eelemCs,
+ eelemNR
+};
+
+enum {
+ emlH, emlF, emlCl, emlBr, emlI,
+ emlCTE, emlCTR, emlCBR, emlCDI,
+ emlNTE, emlNTR2, emlNPI2, emlNDI,
+ emlOTE, emlOTR4, emlOPI2,
+ emlSTE, emlSTR4, emlSPI2,
+ emlPTE, emlSi, emlNR
+};
+
+enum {
+ eqmHF, eqmMP2DZ, eqmMP2TZ, eqmDFT, eqmDFTD,
+ eqmMiller, eqmKang, eqmBosque, eqmSpoel, eqmNR
+};
+#define NQUANT (eqmDFTD+1)
+
+/* Global variable! */
+extern char *lbasis[eqmNR];
+
+#define ecNR 20
+#define eatNR 29
+#define eatExtra 20
+
+#define nullify(x) memset(&(x),0,sizeof(x))
+#define snew(x,n) { (x) = calloc((n),sizeof(*x)); }
+#define srenew(x,n) { if ((x)) { (x) = realloc((x),(n)*sizeof(*x)); } else { snew((x),(n)); } }
+#define sfree(x) { if ((x!=NULL)) { free(x); x = NULL; } }
+
+typedef struct {
+ int nr;
+ char *formula,*molname;
+ int category[ecNR];
+ double weight;
+ int nexperiment;
+ double *experiment,*error;
+ char **pname,**reference;
+ double qm[eqmNR];
+ int frag_comp[eatNR+eatExtra];
+ int elem_comp[eelemNR];
+ int emil_comp[emlNR];
+} t_molprop;
+
+typedef struct {
+ char *name,*atom,*miller,*bosque;
+ int nh,maxbond,hybrid;
+ double train,all,blen;
+} t_spoel;
+
+typedef struct {
+ char *name;
+ int nelec;
+ double value,mass;
+} t_bosque;
+
+typedef struct {
+ char *name;
+ int nelec;
+ double tau_ahc,alpha_ahp;
+} t_miller;
+
+/* Global variables! */
+extern t_miller miller[emlNR];
+extern t_spoel spoel[eatNR+eatExtra];
+extern t_bosque bosque[eelemNR+1];
+extern char *ec_name[ecNR];
+
+extern void write_molprops(char *fn,int npd,t_molprop pd[]);
+
+extern int read_molprops(char *fn,t_molprop **pd,int update_bm);
+/* Return number of molprops. If update_bm != 0 then the Bosque and Miller
+ compositions will be updated based on the Spoel composition */
+
+extern void gaussj(double **a,int n,double **b,int m);
+
+extern void fatal_error(char *str,char *val);
+
+extern int mp_num_prop(t_molprop *mp,char *prop);
+#define mp_num_polar(mp) mp_num_prop(mp,"Polarizability")
+#define mp_num_dip(mp) mp_num_prop(mp,"Dipole")
+
+extern double mp_get_prop(t_molprop *mp,char *prop,int index);
+#define mp_get_polar(mp) mp_get_prop(mp,"Polarizability",0)
+#define mp_get_dip(mp) mp_get_prop(mp,"Dipole",0)
+
+extern char *mp_get_ref_prop(t_molprop *mp,char *prop,int index);
+#define mp_get_ref_polar(mp) mp_get_ref_prop(mp,"Polarizability",0)
+#define mp_get_ref_dip(mp) mp_get_ref_prop(mp,"Dipole",0)
+
+
+extern t_molprop *merge_xml(int argc,char *argv[],char *outf,
+ char *sorted,int *nmolprop);
+
+#endif
static bool bUsePDBcharge = FALSE,bVerbose=FALSE;
static char *molnm = "ICE";
static char *ff = "select";
- static char *qgen[] = { NULL, "None", "Linear", "Yang", "Bultinck", "SM", NULL };
+ static char *qgen[] = { NULL, "None", "Linear", "Yang", "Bultinck",
+ "SM1", "SM2", "SM3", "SM4", NULL };
t_pargs pa[] = {
{ "-ff", FALSE, etSTR, {&ff},
"Select the force field for your simulation." },
alg = eqgNone;
if (qgen[0]) {
- for(alg=1; (alg<=eqgNR); alg++) {
- if (strcmp(qgen[0],qgen[alg]) == 0)
- break;
- }
- if (alg > eqgNR)
+ alg = name2eemtype(qgen[0]);
+ if (alg == -1)
alg = eqgNone;
- else
- alg--;
}
if (alg == eqgNone) {
if (nqa == atoms->nr) {
fp = ftp2FILE(efTOP,NFILE,fnm,"w");
print_top_header(fp,ftp2fn(efTOP,NFILE,fnm),
"Generated by x2top",TRUE, forcefield,1.0);
- /*switch (alg) {
- case eqgYang:
- please_cite(fp,"Yang2006b");
- break;
- case eqgBultinck:
- please_cite(fp,"Bultinck2002a");
- break;
- default:
- break;
- } */
write_top(fp,NULL,mymol.name,atoms,bts,plist,excls,atype,
cgnr,nexcl);
print_top_mols(fp,mymol.name,NULL,0,NULL,1,&mymol);
} t_eemrecord;
static char *eemtype_name[eqgNR] = {
- "None", "Linear", "Yang", "Bultinck", "SM"
+ "None", "Linear", "Yang", "Bultinck", "SM1", "SM2", "SM3", "SM4"
};
-static int name2eemtype(char *name)
+int name2eemtype(char *name)
{
int i;
return -1;
}
-void *read_eemprops(char *fn)
+void *read_eemprops(char *fn,int eemtype)
{
t_eemrecord *eem=NULL;
- char buf[STRLEN],**strings;
- int i,n;
+ char buf[STRLEN],**strings,*ptr;
+ int i,n,nn=0;
char nmbuf[32],algbuf[32];
int elem,row;
double J0,radius,chi0;
snew(eem,1);
snew(eem->eep,n);
for(i=0; (i<n); i++) {
- if (sscanf(strings[i],"%s%s%d%d%lf%lf%lf",nmbuf,algbuf,&elem,&row,
- &J0,&radius,&chi0) != 7)
- gmx_fatal(FARGS,"Error in %s on line %d",buf,i+1);
- eem->eep[i].name = strdup(nmbuf);
- if ((eem->eep[i].eemtype = name2eemtype(algbuf)) == -1)
- gmx_fatal(FARGS,"Error in %s on line %d, unknown algorithm '%s'",
+ ptr = strings[i];
+ while (*ptr && isspace(*ptr))
+ ptr++;
+ if (((ptr) && (*ptr != ';')) &&
+ (sscanf(strings[i],"%s%s%d%d%lf%lf%lf",nmbuf,algbuf,&elem,&row,
+ &J0,&radius,&chi0) == 7)) {
+ if ((eem->eep[nn].eemtype = name2eemtype(algbuf)) == -1)
+ fprintf(stderr,"Warning in %s on line %d, unknown algorithm '%s'\n",
buf,i+1,algbuf);
- eem->eep[i].elem = elem;
- eem->eep[i].row = row;
- eem->eep[i].J0 = J0;
- eem->eep[i].radius = radius;
- eem->eep[i].chi0 = chi0;
+ else if ((eemtype == -1) || (eem->eep[nn].eemtype == eemtype)) {
+ eem->eep[nn].name = strdup(nmbuf);
+ eem->eep[nn].elem = elem;
+ eem->eep[nn].row = row;
+ eem->eep[nn].J0 = J0;
+ eem->eep[nn].radius = radius;
+ eem->eep[nn].chi0 = chi0;
+ nn++;
+ }
+ }
}
- eem->nep = n;
+ eem->nep = nn;
}
+
return eem;
}
t_eemrecord *er = (t_eemrecord *) eem;
int i;
+ fprintf(fp,"; Atom Model Nr Row J_aa w_a Chi_a\n");
for(i=0; (i<er->nep); i++)
fprintf(fp,"%-5s %10s %3d %3d %10.3f %10.3f %10.3f\n",
er->eep[i].name,eemtype_name[er->eep[i].eemtype],
er->eep[i].radius,er->eep[i].chi0);
}
-int eem_getnumprops(void *eem)
+int eem_get_numprops(void *eem)
{
t_eemrecord *er = (t_eemrecord *) eem;
return er->nep;
}
-int eem_getindex(void *eem,char *resname,char *aname,int eemtype)
+int eem_get_index(void *eem,char *resname,char *aname,int eemtype)
{
t_eemrecord *er = (t_eemrecord *) eem;
int i;
else
*wj = 10*(3/(4*er->eep[index].radius));
}
- else if (er->eep[index].eemtype == eqgSM)
+ else if ((er->eep[index].eemtype == eqgSM2) ||
+ (er->eep[index].eemtype == eqgSM3) ||
+ (er->eep[index].eemtype == eqgSM4))
*wj = 10.0/er->eep[index].radius;
else
*wj = 0;
real eem_get_j00(void *eem,char *resname,char *aname,real *wj,real qH,int eemtype)
{
- int k = eem_getindex(eem,resname,aname,eemtype);
+ int k = eem_get_index(eem,resname,aname,eemtype);
return lo_get_j00(eem,k,wj,qH);
}
return er->eep[index].radius;
}
+real eem_get_elem(void *eem,int index)
+{
+ t_eemrecord *er = (t_eemrecord *) eem;
+
+ range_check(index,0,er->nep);
+
+ return er->eep[index].elem;
+}
+
void eem_set_props(void *eem,int index,real J0,real radius,real chi0)
{
t_eemrecord *er = (t_eemrecord *) eem;
sfree(a);
}
-void mat_mult(FILE *fp,int n,int m,double **x,double **y,double **z)
+void matrix_multiply(FILE *fp,int n,int m,double **x,double **y,double **z)
{
int i,j,k;
}
}
-void mat_inv(FILE *fp,int n,double **a)
+void matrix_invert(FILE *fp,int n,double **a)
{
int i,j,m,lda,*ipiv,lwork,info;
double **test,**id,*work;
}
if (fp) {
id = alloc_matrix(n,n);
- mat_mult(fp,n,n,test,a,id);
+ matrix_multiply(fp,n,n,test,a,id);
fprintf(fp,"And here is the product of A and Ainv\n");
for(i=0; (i<n); i++) {
for(j=0; (j<n); j++)
typedef struct {
int natom,eemtype;
int *index; /* In the Yang array */
- real *chi,*chi0,*qq,*wj,qtotal;
+ int *elem;
+ real *chi,*chi0,*rhs,*qq,*wj,qtotal;
real **Jab;
rvec *x;
- real chiav;
+ real chieq;
} t_qgen;
static real coul_slater_slater(real w,real r)
return 1/r;
}
-static real calc_jab(rvec xi,rvec xj,real wi,real wj)
+static real calc_jab(rvec xi,rvec xj,real wi,real wj,int eemtype)
{
rvec dx;
real r,wij;
- real e0=0,e1=0,e2=0;
+ real eNN=0,eSS=0;
rvec_sub(xi,xj,dx);
r = norm(dx);
+ if (r == 0)
+ gmx_fatal(FARGS,"Zero distance between atoms!\n");
+
+ switch (eemtype) {
+ case eqgBultinck:
+ case eqgSM1:
+ eNN = coul_nucl_nucl(0,r);
+ break;
+ case eqgSM2:
+ case eqgSM3:
+ case eqgSM4:
+ if ((wi > 0) && (wj > 0)) {
+ wij = 2*(wi + wj)/(wi*wj);
+ eSS = coul_slater_slater(wij,r);
+ }
+ else {
+ eNN = coul_nucl_nucl(0,r);
+ }
+ break;
+ case eqgYang:
+ default:
+ gmx_fatal(FARGS,"Can not treat algorithm %d yet in calc_jab",eemtype);
+ }
+
+ return ONE_4PI_EPS0*(eNN+eSS)/ELECTRONVOLT;
+}
+
+static real calc_j1(rvec xi,rvec xj,real wi,real wj,int eemtype)
+{
+ rvec dx;
+ real r,wij;
+ real eNN=0,eNS=0,eSN=0,eSS=0;
+
+ rvec_sub(xi,xj,dx);
+ r = norm(dx);
+ if (r == 0)
+ gmx_fatal(FARGS,"Zero distance between atoms!\n");
- e0 = coul_nucl_nucl(0,r);
if ((wi > 0) && (wj > 0)) {
- wij = 2/(1/wi + 1/wj); /* Dit kan geoptimaliseerd worden */
- e1 = coul_slater_nucl(wij,r);
- e2 = coul_slater_slater(wij,r);
+ wij = 2*(wi + wj)/(wi*wj);
+ eSN = coul_slater_nucl(wj,r);
+ eSS = coul_slater_slater(wij,r);
}
-
- return ONE_4PI_EPS0*(e0-2*e1+e2)/ELECTRONVOLT;
+ else {
+ eNN = coul_nucl_nucl(0,r);
+ }
+ return ONE_4PI_EPS0*(eNN+eSN-eSS)/ELECTRONVOLT;
}
static real get_chi0(void *atomprop,char *resnm,char *name)
static void solve_q_eem(FILE *fp,t_qgen *qgen,real hardness_factor)
{
- double **a,**b,qtot,chieq;
+ double **a,**b,qtot;
int i,j,n,nn;
n = qgen->natom+1;
for(j=0; (j<n); j++)
b[i][j] = a[i][j];
}
- mat_inv(fp,n,a);
+ matrix_invert(fp,n,a);
qtot = 0;
for(i=0; (i<n-1); i++) {
qgen->qq[i] = 0;
for(j=0; (j<n-1); j++) {
- qgen->qq[i] += -a[i][j]*qgen->chi0[j];
+ qgen->qq[i] += a[i][j]*qgen->rhs[j];
}
qtot += qgen->qq[i];
}
- if (fp) {
- chieq = 0;
- for(i=0; (i<n-1); i++) {
- qgen->chi[i] = 0;
- for(j=0; (j<n-1); j++) {
- qgen->chi[i] += b[i][j]*qgen->qq[j];
- }
- chieq += qgen->chi[i];
- }
- qgen->chiav = chieq/qgen->natom;
- free_matrix(b,n);
- }
+ qgen->chieq = 0;
+ for(j=0; (j<n); j++)
+ qgen->chieq += a[n-1][j]*qgen->rhs[j];
- if (fabs(qtot - qgen->qtotal) > 1e-3)
+ if (fabs(qtot - qgen->qtotal) > 1e-2)
fprintf(stderr,"qtot = %g, it should be %g\n",qtot,qgen->qtotal);
free_matrix(a,n);
}
-static void qgen_calc_Jab(t_qgen *qgen,void *eem)
+static void qgen_calc_Jab(t_qgen *qgen,void *eem,int eemtype)
{
int i,j;
double wi,wj;
qgen->Jab[i][i] = lo_get_j00(eem,qgen->index[i],&(qgen->wj[i]),qgen->qq[i]);
}
for(i=0; (i<qgen->natom); i++) {
- wi = qgen->wj[i];
+ wi = qgen->wj[i];
+ qgen->rhs[i] = -qgen->chi0[i];
for(j=0; (j<qgen->natom); j++) {
if (i != j) {
wj = qgen->wj[j];
- qgen->Jab[i][j] = qgen->Jab[j][i] =
- calc_jab(qgen->x[i],qgen->x[j],wi,wj);
+ qgen->Jab[i][j] = calc_jab(qgen->x[i],qgen->x[j],wi,wj,qgen->eemtype);
+ if ((eemtype == eqgSM3) || (eemtype == eqgSM4))
+ qgen->rhs[i] -= qgen->elem[j]*calc_j1(qgen->x[i],qgen->x[j],wi,wj,qgen->eemtype);
}
}
}
qgen->natom = atoms->nr;
qgen->eemtype = eemtype;
snew(qgen->chi0,atoms->nr);
+ snew(qgen->rhs,atoms->nr);
+ snew(qgen->elem,atoms->nr);
snew(qgen->chi,atoms->nr);
snew(qgen->Jab,atoms->nr);
snew(qgen->wj,atoms->nr);
qgen->x = x;
for(i=0; (i<atoms->nr); i++) {
snew(qgen->Jab[i],atoms->nr);
- qgen->index[i] = eem_getindex(eem,
- *(atoms->resname[atoms->atom[i].resnr]),
- *(atoms->atomname[i]),qgen->eemtype);
+ qgen->index[i] = eem_get_index(eem,
+ *(atoms->resname[atoms->atom[i].resnr]),
+ *(atoms->atomname[i]),qgen->eemtype);
+ qgen->elem[i] = eem_get_elem(eem,qgen->index[i]);
if (qgen->index[i] == -1)
- gmx_fatal(FARGS,"Can not find index for %s %s",
+ gmx_fatal(FARGS,"Can not find index for %s %s. Eemtype = %d",
*(atoms->resname[atoms->atom[i].resnr]),
- *(atoms->atomname[i]));
+ *(atoms->atomname[i]),eemtype);
qgen->chi0[i] = eem_get_chi0(eem,qgen->index[i]);
}
}
}
if (fp)
- fprintf(fp,"<chieq> = %10g\n",qgen->chiav);
+ fprintf(fp,"<chieq> = %10g\n",qgen->chieq);
sfree(qgen->chi0);
sfree(qgen->chi);
sfree(qgen->wj);
qq[i] = qgen->qq[i];
iter = 0;
do {
- qgen_calc_Jab(qgen,eem);
+ qgen_calc_Jab(qgen,eem,eqgYang);
solve_q_eem(debug,qgen,2.0);
rms = 0;
for(i=0; (i<atoms->nr); i++) {
done_qgen(stdout,atoms,qgen);
}
-void generate_charges_sm(FILE *fp,
+real generate_charges_sm(FILE *fp,
void *eem,t_atoms *atoms,rvec x[],
real tol,int maxiter,void *atomprop,
- real qtotref)
+ real qtotref,int eemtype)
{
/* Use Rappe and Goddard derivative for now */
t_qgen *qgen;
- real *qq;
+ real *qq,chieq;
int i,iter;
real rms,mu;
if (fp)
fprintf(fp,"Generating charges using Van der Spoel & Van Maaren algorithm\n");
- qgen = init_qgen(eem,atoms,atomprop,x,eqgSM);
+ qgen = init_qgen(eem,atoms,atomprop,x,eemtype);
snew(qq,atoms->nr);
for(i=0; (i<atoms->nr); i++)
qq[i] = qgen->qq[i];
iter = 0;
do {
- qgen_calc_Jab(qgen,eem);
+ qgen_calc_Jab(qgen,eem,eemtype);
solve_q_eem(fp,qgen,2.0);
rms = 0;
for(i=0; (i<atoms->nr); i++) {
fprintf(fp,"Did not converge with %d iterations. RMS = %g\n",maxiter,rms);
}
sfree(qq);
+ chieq = qgen->chieq;
done_qgen(fp,atoms,qgen);
sfree(qgen);
+
+ return chieq;
}
static void generate_charges_bultinck(void *eem,t_atoms *atoms,rvec x[],
printf("Generating charges using Bultinck algorithm\n");
qgen = init_qgen(eem,atoms,atomprop,x,eqgBultinck);
- qgen_calc_Jab(qgen,eem);
+ qgen_calc_Jab(qgen,eem,eqgBultinck);
solve_q_eem(NULL,qgen,2.0);
done_qgen(stdout,atoms,qgen);
sfree(chi0);
}
-void assign_charge_alpha(int alg,t_atoms *atoms,rvec x[],
+void assign_charge_alpha(int eemtype,t_atoms *atoms,rvec x[],
t_params *bonds,real tol,real fac,int maxiter,
void *atomprop,real qtotref)
{
int i;
void *eem;
- eem = read_eemprops(NULL);
+ eem = read_eemprops(NULL,-1);
if (debug)
write_eemprops(debug,eem);
- if ((eem == NULL) && (alg > eqgLinear))
+ if ((eem == NULL) && (eemtype > eqgLinear))
gmx_fatal(FARGS,"Nothing interesting in eemprops.dat");
/* Generate charges */
- switch (alg) {
+ switch (eemtype) {
case eqgNone:
for(i=0; (i<atoms->nr); i++) {
atoms->atom[i].q = atoms->atom[i].qB = 0;
please_cite(stdout,"Bultinck2002a");
generate_charges_bultinck(eem,atoms,x,tol,maxiter,atomprop);
break;
- case eqgSM:
- generate_charges_sm(stdout,eem,atoms,x,tol,maxiter,atomprop,qtotref);
+ case eqgSM1:
+ case eqgSM2:
+ case eqgSM3:
+ case eqgSM4:
+ (void) generate_charges_sm(stdout,eem,atoms,x,tol,maxiter,atomprop,
+ qtotref,eemtype);
break;
default:
- gmx_fatal(FARGS,"Algorithm %d out of range in assign_charge_alpha",alg);
+ gmx_fatal(FARGS,"Algorithm %d out of range in assign_charge_alpha",
+ eemtype);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff