src/gmxlib/xmlio.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #include <stdlib.h>
  37 #include "typedefs.h"
  38 #ifdef HAVE_XML
  39 #include <libxml/parser.h>
  40 #include <libxml/tree.h>
  41 #include "fatal.h"
  42 #include "string.h"
  43 #include "futil.h"
  44 #include "smalloc.h"
  45 #include "names.h"
  46 #include "assert.h"
  47 #include "symtab.h"
  48 #include "macros.h"
  49 #include "symtab.h"
  50 #include "xmlio.h"
  51
  52 static const char *xyz_names[]    = {
  53   "x", "y", "z"
  54 };
  55 static const char *tensor_names[] = {
  56   "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"
  57 };
  58
  59 typedef struct {
  60   int        nx,nv,nf,ntop,nbox,ninprec;
  61   int        step,natoms;
  62   real       t,lambda;
  63   t_inputrec *ir;
  64   rvec       *box;
  65   rvec       *x,*v,*f;
  66   t_topology *top;
  67 } t_xmlrec;
  68
  69 typedef struct {
  70   char *name;
  71   real value;
  72 } t_masstype;
  73
  74 static const char *xmltypes[] = {
  75   NULL,
  76   "XML_ELEMENT_NODE",
  77   "XML_ATTRIBUTE_NODE",
  78   "XML_TEXT_NODE",
  79   "XML_CDATA_SECTION_NODE",
  80   "XML_ENTITY_REF_NODE",
  81   "XML_ENTITY_NODE",
  82   "XML_PI_NODE",
  83   "XML_COMMENT_NODE",
  84   "XML_DOCUMENT_NODE",
  85   "XML_DOCUMENT_TYPE_NODE",
  86   "XML_DOCUMENT_FRAG_NODE",
  87   "XML_NOTATION_NODE",
  88   "XML_HTML_DOCUMENT_NODE",
  89   "XML_DTD_NODE",
  90   "XML_ELEMENT_DECL",
  91   "XML_ATTRIBUTE_DECL",
  92   "XML_ENTITY_DECL",
  93   "XML_NAMESPACE_DECL",
  94   "XML_XINCLUDE_START",
  95   "XML_XINCLUDE_END"
  96 };
  97 #define NXMLTYPES asize(xmltypes)
  98
  99 extern int xmlDoValidityCheckingDefaultValue;
 100
 101 enum {
 102   exmlGROMACS,
 103   exmlINPUTREC,    exmlOUTPUT,      exmlCOUPLING,    exmlCUTOFF,
 104   exmlPMEPARM,     exmlTCOUPLING,   exmlPCOUPLING,   exmlTCPARM,
 105   exmlREFPRES,     exmlCOMPRESS,    exmlRVEC,
 106   exmlSYSTEM,      exmlCOMPOSITION, exmlMOLECULE,    exmlATOM,
 107   exmlTOPOLOGY,    exmlBONDS,       exmlANGLES,      exmlDIHEDRALS,
 108   exmlFORCEFIELD,  exmlMASSTYPE,    exmlBOX,
 109   exmlCOORDINATES, exmlVELOCITIES,  exmlFORCES,      exmlNR
 110 };
 111
 112 static const char *exml_names[] = {
 113   "gromacs",
 114   /* Inputrec stuff */
 115   "parameters",  "output",     "coupling", "cutoff",   "pmeparm",
 116   "tcoupling",   "pcoupling",  "tcparm",   "p-ref",    "compressibility", "rvec",
 117   /* System description */
 118   "system",      "composition", "molecule","atom",
 119   /* Topology description */
 120   "topology",
 121   "bonds",       "angles",      "dihedrals",
 122   /* Force field */
 123   "forcefield",  "masstype",    "cell",
 124   /* Coordinates etc. */
 125   "coordinates", "velocities", "forces"
 126 };
 127
 128 static int find_elem(char *name,int nr,const char *names[])
 129 {
 130   int i;
 131
 132   for(i=0; (i<nr); i++)
 133     if (strcmp(name,names[i]) == 0)
 134       break;
 135   if (i == nr)
 136     fatal_error(0,"Unknown element name %s",name);
 137
 138   return i;
 139 }
 140
 141 static char *sp(int n, char buf[], int maxindent)
 142 {
 143   int i;
 144   if(n>=maxindent)
 145     n=maxindent-1;
 146
 147   /* Don't indent more than maxindent characters */
 148   for(i=0; (i<n); i++)
 149     buf[i] = ' ';
 150   buf[i] = '\0';
 151
 152   return buf;
 153 }
 154
 155 static void process_attr(FILE *fp,xmlAttrPtr attr,int elem,
 156                          int indent,t_xmlrec *xml)
 157 {
 158   char *attrname,*attrval;
 159   char buf[100];
 160
 161   while (attr != NULL) {
 162     attrname = (char *)attr->name;
 163     attrval  = (char *)attr->children->content;
 164
 165 #define atest(s) ((strcmp(attrname,s) == 0) && (attrval != NULL))
 166     switch (elem) {
 167     case exmlGROMACS:
 168       if (atest("title"))
 169         xml->top->name = put_symtab(&xml->top->symtab,attrval);
 170       break;
 171     case exmlINPUTREC:
 172       if (atest("algorithm"))
 173         xml->ir->eI = find_elem(attrval,eiNR,ei_names);
 174       break;
 175     case exmlOUTPUT:
 176     case exmlCOUPLING:
 177     case exmlCUTOFF:
 178     case exmlPMEPARM:
 179     case exmlTCOUPLING:
 180     case exmlPCOUPLING:
 181     case exmlTCPARM:
 182     case exmlREFPRES:
 183     case exmlCOMPRESS:
 184     case exmlRVEC:
 185     case exmlSYSTEM:
 186     case exmlCOMPOSITION:
 187     case exmlMOLECULE:
 188     case exmlATOM:
 189     case exmlTOPOLOGY:
 190     case exmlBONDS:
 191     case exmlANGLES:
 192     case exmlDIHEDRALS:
 193     case exmlFORCEFIELD:
 194     case exmlMASSTYPE:
 195     case exmlCOORDINATES:
 196     case exmlVELOCITIES:
 197     case exmlFORCES:
 198     default:
 199       if (fp)
 200         fprintf(fp,"%sProperty: '%s' Value: '%s'\n",sp(indent,buf,99),
 201                 attrname,attrval);
 202     }
 203     attr = attr->next;
 204 #undef atest
 205   }
 206 }
 207
 208 static void process_tree(FILE *fp,xmlNodePtr tree,int indent,t_xmlrec *xml)
 209 {
 210   int elem;
 211   char buf[100];
 212
 213   while (tree != NULL) {
 214     switch (tree->type) {
 215     case XML_ELEMENT_NODE:
 216       elem = find_elem((char *)tree->name,exmlNR,exml_names);
 217       if (fp)
 218         fprintf(fp,"%sElement node name %s\n",sp(indent,buf,99),(char *)tree->name);
 219
 220       process_attr(fp,tree->properties,elem,indent+2,xml);
 221
 222       if (tree->children)
 223         process_tree(fp,tree->children,indent+2,xml);
 224       break;
 225     case XML_COMMENT_NODE:
 226       if (fp)
 227         fprintf(fp,"Comment node encountered\n");
 228       break;
 229     case XML_PI_NODE:
 230     case XML_TEXT_NODE:
 231       /* Silently ignore these for now */
 232       break;
 233     default:
 234       if (fp) {
 235         if ((tree->type > 0) && (tree->type < NXMLTYPES))
 236           fprintf(fp,"Node type %s encountered with name %s\n",
 237                   xmltypes[tree->type],(char *)tree->name);
 238         else
 239           fprintf(fp,"Node type %d encountered\n",tree->type);
 240       }
 241     }
 242     tree = tree->next;
 243   }
 244 }
 245
 246 void read_xml(char *fn,int *step,real *t,real *lambda,
 247               t_inputrec *ir,rvec *box,int *natoms,
 248               rvec **x,rvec **v,rvec **f,t_topology *top)
 249 {
 250   xmlDocPtr  doc;
 251   t_xmlrec   *xml;
 252
 253   xmlDoValidityCheckingDefaultValue = 1;
 254   assert(asize(exml_names) == exmlNR);
 255   if ((doc = xmlParseFile(fn)) == NULL)
 256     fatal_error(0,"Reading XML file %s. Run a syntax checker such as nsgmls.",
 257                 fn);
 258
 259   snew(xml,1);
 260   xml->ir     = ir;
 261   xml->box    = box;
 262   xml->top    = top;
 263   process_tree(debug,doc->children,0,xml);
 264
 265   xmlFreeDoc(doc);
 266   sfree(xml);
 267 }
 268
 269 static void add_xml_int(xmlNodePtr ptr,const char *name,int val)
 270 {
 271   char buf[32];
 272   sprintf(buf,"%d",val);
 273   if (xmlSetProp(ptr,name,buf) == 0)
 274     fatal_error(0,"Setting %s %d",name,val);
 275 }
 276
 277 static void add_xml_real(xmlNodePtr ptr,const char *name,real val)
 278 {
 279   char buf[32];
 280   sprintf(buf,"%g",val);
 281   if (xmlSetProp(ptr,name,buf) == 0)
 282     fatal_error(0,"Setting %s %f",name,val);
 283 }
 284
 285 static void add_xml_rvec(xmlNodePtr parent,int id,rvec val)
 286 {
 287   xmlNodePtr rvptr;
 288   int  m;
 289   char buf[32];
 290
 291   if ((rvptr = xmlNewChild(parent,NULL,"rvec",NULL)) == NULL)
 292     fatal_error(0,"Creating rvec element");
 293   add_xml_int(rvptr,"id",id);
 294   for(m=0; (m<DIM); m++) {
 295     sprintf(buf,"%g",val[m]);
 296     if (xmlSetProp(rvptr,xyz_names[m],buf) == 0)
 297       fatal_error(0,"Setting %s %f",xyz_names[m],val[m]);
 298   }
 299 }
 300
 301 static void add_xml_tensor(xmlNodePtr parent,tensor val)
 302 {
 303   xmlNodePtr tptr;
 304   int  m,n;
 305
 306   if ((tptr = xmlNewChild(parent,NULL,"tensor",NULL)) == NULL)
 307     fatal_error(0,"Creating tensor element");
 308
 309   for(m=0; (m<DIM); m++)
 310     add_xml_real(tptr,tensor_names[m*DIM+m],val[m][m]);
 311   if ((val[XX][YY] != 0) || (val[XX][ZZ] != 0) ||
 312       (val[YY][XX] != 0) || (val[YY][ZZ] != 0) ||
 313       (val[ZZ][XX] != 0) || (val[ZZ][YY] != 0)) {
 314     for(m=0; (m<DIM); m++)
 315       for(n=0; (n<DIM); n++)
 316         if (m != n)
 317           add_xml_real(tptr,tensor_names[m*DIM+n],val[m][n]);
 318   }
 319 }
 320
 321 static void add_xml_char(xmlNodePtr ptr,const char *name,const char *val)
 322 {
 323   if (xmlSetProp(ptr,name,val) == 0)
 324     fatal_error(0,"Setting %s %s",name,val);
 325 }
 326
 327 static xmlNodePtr add_xml_child(xmlNodePtr parent,int type)
 328 {
 329   xmlNodePtr child;
 330
 331   if ((child = xmlNewChild(parent,NULL,exml_names[type],NULL)) == NULL)
 332     fatal_error(0,"Creating %s element",exml_names[type]);
 333
 334   return child;
 335 }
 336
 337 static xmlNodePtr add_xml_comment(xmlDocPtr doc,
 338                                   xmlNodePtr prev,char *comment)
 339 {
 340   xmlNodePtr comm,ptr;
 341
 342   if ((comm = xmlNewComment((xmlChar *)comment)) == NULL)
 343     fatal_error(0,"Creating doc comment element");
 344   ptr = prev;
 345   while (ptr->next != NULL)
 346     ptr=ptr->next;
 347   ptr->next    = comm;
 348   comm->prev   = ptr;
 349   comm->doc    = doc;
 350
 351   return comm;
 352 }
 353
 354 static void add_xml_inputrec(xmlNodePtr parent,t_inputrec *ir,t_atoms *atoms)
 355 {
 356   int        i;
 357   xmlNodePtr irptr,outputptr,tcptr,tcparm,pcptr,refpres,compress;
 358   xmlNodePtr cutoffptr,pmeptr;
 359
 360   irptr = add_xml_child(parent,exmlINPUTREC);
 361   add_xml_char(irptr,"algorithm",ei_names[ir->eI]);
 362   add_xml_int(irptr,"nsteps",ir->nsteps);
 363   add_xml_real(irptr,"init-t",ir->init_t);
 364   add_xml_real(irptr,"delta-t",ir->delta_t);
 365
 366   outputptr = add_xml_child(irptr,exmlOUTPUT);
 367   add_xml_int(outputptr,"log",ir->nstlog);
 368   add_xml_int(outputptr,"x-trr",ir->nstxout);
 369   add_xml_int(outputptr,"v-trr",ir->nstvout);
 370   add_xml_int(outputptr,"f-trr",ir->nstfout);
 371   add_xml_int(outputptr,"energy",ir->nstenergy);
 372   add_xml_int(outputptr,"x-xtc",ir->nstxtcout);
 373   add_xml_int(outputptr,"xtc-precision",ir->xtcprec);
 374   add_xml_int(outputptr,"andersen_seed",ir->andersen_seed);
 375
 376   tcptr = add_xml_child(irptr,exmlTCOUPLING);
 377   add_xml_char(tcptr,"algorithm",etcoupl_names[ir->etc]);
 378
 379   assert(ir->opts.ngtc == atoms->grps[egcTC].nr);
 380   for(i=0; (i<ir->opts.ngtc); i++) {
 381     tcparm = add_xml_child(tcptr,exmlTCPARM);
 382     add_xml_char(tcparm,"groupname",*atoms->grpname[atoms->grps[egcTC].nm_ind[i]]);
 383     add_xml_real(tcparm,"t-ref",ir->opts.ref_t[i]);
 384     add_xml_real(tcparm,"tau-t",ir->opts.tau_t[i]);
 385   }
 386
 387   pcptr = add_xml_child(irptr,exmlPCOUPLING);
 388   add_xml_char(pcptr,"algorithm",epcoupl_names[ir->epc]);
 389   add_xml_char(pcptr,"type",epcoupltype_names[ir->epct]);
 390   add_xml_real(pcptr,"tau-p",ir->tau_p);
 391
 392   refpres = add_xml_child(pcptr,exmlREFPRES);
 393   add_xml_tensor(refpres,ir->ref_p);
 394
 395   compress = add_xml_child(pcptr,exmlCOMPRESS);
 396   add_xml_tensor(compress,ir->compress);
 397
 398   cutoffptr = add_xml_child(irptr,exmlCUTOFF);
 399   add_xml_real(cutoffptr,"rlist",ir->rlist);
 400   add_xml_real(cutoffptr,"rvdw",ir->rvdw);
 401   add_xml_real(cutoffptr,"rcoulomb",ir->rcoulomb);
 402   add_xml_real(cutoffptr,"rcoulswitch",ir->rcoulomb_switch);
 403   add_xml_real(cutoffptr,"rvdwswitch",ir->rvdw_switch);
 404   add_xml_real(cutoffptr,"epsilonr",ir->epsilon_r);
 405   add_xml_int(cutoffptr,"nstlist",ir->nstlist);
 406   add_xml_char(cutoffptr,"nstype",ens_names[ir->ns_type]);
 407   add_xml_char(cutoffptr,"domdecomp",yesno_names[ir->bDomDecomp]);
 408   add_xml_char(cutoffptr,"decompdir",xyz_names[ir->decomp_dir]);
 409   add_xml_char(cutoffptr,"coulombtype",eel_names[ir->coulombtype]);
 410   add_xml_char(cutoffptr,"vdwtype",evdw_names[ir->vdwtype]);
 411   if (ir->coulombtype == eelPME) {
 412     pmeptr = add_xml_child(cutoffptr,exmlPMEPARM);
 413     add_xml_int(pmeptr,"nkx",ir->nkx);
 414     add_xml_int(pmeptr,"nky",ir->nky);
 415     add_xml_int(pmeptr,"nkz",ir->nkz);
 416     add_xml_int(pmeptr,"pmeorder",ir->pme_order);
 417     add_xml_real(pmeptr,"ewaldrtol",ir->ewald_rtol);
 418     add_xml_real(pmeptr,"epssurface",ir->epsilon_surface);
 419     add_xml_char(pmeptr,"optfft",yesno_names[ir->bOptFFT]);
 420   }
 421 }
 422
 423 static void add_xml_molecule(xmlNodePtr parent,t_atoms *atoms,
 424                              int nmt,t_masstype mt[])
 425 {
 426   xmlNodePtr ptr;
 427
 428   ptr = add_xml_child(parent,exmlMOLECULE);
 429 }
 430
 431 static void add_xml_idef(xmlNodePtr parent,t_idef *idef)
 432 {
 433
 434 }
 435
 436 static void add_xml_rvecs(xmlNodePtr parent,int type,int natoms,rvec *xvf)
 437 {
 438   xmlNodePtr xptr;
 439   int i;
 440
 441   xptr = add_xml_child(parent,type);
 442   for(i=0; (i<natoms); i++)
 443     add_xml_rvec(xptr,i+1,xvf[i]);
 444 }
 445
 446 static t_masstype *mk_masstype(int nmol,t_atoms atoms[],int *nmt)
 447 {
 448   int        i,j,k,nm;
 449   t_masstype *mt=NULL;
 450
 451   nm = 0;
 452   for(i=0; (i<nmol); i++) {
 453     for(j=0; (j<atoms[i].nr); j++) {
 454       for(k=0; (k<nm); k++)
 455         if (strcmp(*atoms[i].atomname[j],mt[k].name) == 0)
 456           break;
 457       if (k == nm) {
 458         srenew(mt,nm+1);
 459         mt[nm].name  = strdup(*atoms[i].atomname[j]);
 460         mt[nm].value = atoms[i].atom[j].m;
 461         nm++;
 462       }
 463     }
 464   }
 465   *nmt = nm;
 466
 467   return mt;
 468 }
 469
 470 static void add_xml_masstype(xmlNodePtr parent,int nmt,t_masstype mt[])
 471 {
 472   int i;
 473   xmlNodePtr ptr;
 474
 475   for(i=0; (i<nmt); i++) {
 476     ptr = add_xml_child(parent,exmlMASSTYPE);
 477     add_xml_char(ptr,"name",mt[i].name);
 478     add_xml_real(ptr,"value",mt[i].value);
 479   }
 480 }
 481
 482 void write_xml(char *fn,char *title,t_inputrec *ir,rvec *box,
 483                int natoms,rvec *x,rvec *v,rvec *f,
 484                int nmol,t_atoms atoms[],t_idef *idef)
 485 {
 486   xmlDocPtr  doc;
 487   xmlDtdPtr  dtd;
 488   xmlNodePtr myroot;
 489   t_masstype *mt;
 490   int        i,nmt;
 491   char       *libdtdname,*dtdname;
 492
 493   dtdname    = "gromacs.dtd";
 494   libdtdname = libfn(dtdname);
 495
 496   if ((doc = xmlNewDoc("1.0")) == NULL)
 497     fatal_error(0,"Creating XML document");
 498
 499   if ((dtd = xmlCreateIntSubset(doc,"gromacs",libdtdname,dtdname)) == NULL)
 500     fatal_error(0,"Creating XML DTD");
 501
 502   if ((myroot = xmlNewDocNode(doc,NULL,"gromacs",NULL)) == NULL)
 503     fatal_error(0,"Creating root element");
 504   dtd->next = myroot;
 505   myroot->prev = (xmlNodePtr) dtd;
 506
 507   /* Title of the system */
 508   if (title)
 509     add_xml_char(myroot,"title",title);
 510
 511   /* Add inputrec */
 512   if (ir) {
 513     xmlNodePtr ptr = myroot;
 514     while(ptr->next != NULL)
 515       ptr = ptr->next;
 516     ptr->next = xmlNewDocComment(doc,(xmlChar *)"Here starts parameter section");
 517     ptr->next->prev = ptr;
 518     fprintf(stderr,"Comment type is %s\n",xmltypes[ptr->next->type]);
 519
 520     add_xml_inputrec(myroot,ir,&atoms[0]);
 521   }
 522   /* Generate masstypes */
 523   mt = mk_masstype(nmol,atoms,&nmt);
 524   add_xml_masstype(myroot,nmt,mt);
 525
 526   /* Add molecule definitions */
 527   for(i=0; (i<nmol); i++)
 528     add_xml_molecule(myroot,&(atoms[i]),nmt,mt);
 529
 530   /* Add force field */
 531   if (idef)
 532     add_xml_idef(myroot,idef);
 533
 534   /* Box */
 535   if (box)
 536     add_xml_tensor(add_xml_child(myroot,exmlBOX),box);
 537
 538   /* Coordinates */
 539   if (x)
 540     add_xml_rvecs(myroot,exmlCOORDINATES,natoms,x);
 541
 542   /* Velocities */
 543   if (v)
 544     add_xml_rvecs(myroot,exmlVELOCITIES,natoms,v);
 545
 546   /* Forces */
 547   if (f)
 548     add_xml_rvecs(myroot,exmlFORCES,natoms,f);
 549
 550   xmlSetDocCompressMode(doc,0);
 551   xmlIndentTreeOutput = 1;
 552   if (xmlSaveFormatFileEnc(fn,doc,"ISO-8859-1",2) == 0)
 553     fatal_error(0,"Saving file %s",fn);
 554   xmlFreeDoc(doc);
 555 }
 556
 557 #endif