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