My Documents/Codes/MD++.svn/trunk/src/md.h

Go to the documentation of this file.

/*
  md.h
  by Wei Cai  caiwei@stanford.edu
  Last Modified : Mon Sep  8 17:10:40 2008

  FUNCTION  :  Easy-to-use MD simulation Framefork

  Featuring :  1. Scripting input
               2. X-window display
               3. automatic simulation log
               4. convenient configuration and output file handle
               5. perfect lattice and dislocation config creator
               6. conjugate-gradient relaxation
               7. NVT MD simulation
               8. NPH and NPT simulation (Parrinello-Rahman + Nose-Hoover)
               9. Cell-Verlist combined neighbor list
               
  This is a single CPU code.
  Previous shared-memory option removed.
  May implment MPI parallelization in the future.
*/


#ifndef _MD_H
#define _MD_H

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#ifdef _USETCL
#include "tcl.h"
#ifdef _USETK
#include "tk.h"
#endif
#endif

#ifdef _USEFFTW
#include "fftw3.h"
#else
#define fftw_plan int
#endif

#include "general.h"
#include "organizer.h"
#include "display.h"
#include "linalg3.h"
#include "relax.h"

/* Internal units
   Energy:  eV
   Length:  A
   Stress:  eV/A^3

   Stored Internal Variables
   Position:     r               (in A)
   Velocity:     v*TIMESTEP      (in A)
   Accelaration: a*TIMESTEP^2/2  (in A)
*/

/* Physical constants */
#define AVO  6.0221367E+23
#define EV   1.6021772E-19     /* (J)    */
#define BOLZ 1.380658E-23      /* (J/K) BOLZ=KB*EV */
#define KB   0.8617336E-4      /* (eV/K) */
#define PLANCK_H 4.135667e-15  /* eV/s */

/* mass unit conversion from (gram/mol) to (eV/(A/ps)^2) */
#define MASSCONVERT   1e-3/AVO/EV/1e20*1e24

/* physical constants for Ewald summation */
#define P_CLMB_E 1.0 // for the correction of Coulomb potential (dielectric function constant) 
#define P_CLMB (P_E*P_E/P_E*1.0e10)/(4.0*M_PI*P_EPSILON0*P_CLMB_E) // (eV*A) conversion factor for the Coulomb force
#define P_SQRT_CLMB sqrt(P_CLMB) //3.79470135954879
#define M_SQRT_PI sqrt(M_PI)     //1.77245385090552


#ifndef max
#define max(a,b) (((a)>(b))?(a):(b))
#endif
#ifndef min
#define min(a,b) (((a)<(b))?(a):(b))
#endif

#ifndef SQR
#define SQR(A) ((A)*(A))
#endif

#ifndef CUBE
#define CUBE(A) ((A)*(A)*(A))
#endif

#ifndef POW6
#define POW6(A) SQR(CUBE(A))
#endif

#ifndef POW12
#define POW12(A) SQR(POW6(A))
#endif

#define PME_int_mod(x,y) (((x%y)+y)%y) /* Particle Mesh Ewald: mod function */

#define MAXSPECIES 10          /* maximum number of species */
#define MAXCOLORS  10          /* maximum number of colors in display */
#define MAXCONSTRAINATOMS 2001 /* maximum number of atoms in constrained minimization */
#define MAXCMDLEN  1000        /* maximum length of command string */
#define MAXNHCLEN  20          /* maximum length of NH chain */

/* Random number generator for platforms not supporting drand48 (e.g. cygwin) */
#ifdef _NODRAND48
#define drand48 (1.0/RAND_MAX)*rand
#endif

//shared memory MP deprecated replaced by MPI
//#define MP_Sync() {}
//#define MP_BeginMasterOnly() {}
//#define MP_EndMasterOnly() {}
//#define MP_Free(s) free(s)
//#define MP_cpu() 0
//#define MP_ncpu() 1

#define Realloc(p,t,s) {p=(t *)realloc(p,sizeof(t)*s);}
#define bindcommand_sync bindcommand



/* Input/output files */
class CNFile : public AUXFile /* configuration(cn) file */
{
public:
    CNFile(int t):AUXFile(t){};
    virtual char * describe();
private:
    virtual int writeblock(void *p);
    virtual int readblock(void *p);
};

class POSCARFile : public AUXFile /* POSCAR (vasp) file */
{
public:
    POSCARFile(int t):AUXFile(t){};
    virtual char * describe();
private:
    virtual int writeblock(void *p);
    virtual int readblock(void *p);
};

class LAMMPSFile : public AUXFile /* LAMMPS file */
{
public:
    LAMMPSFile(int t):AUXFile(t){};
    virtual char * describe();
private:
    virtual int writeblock(void *p);
    virtual int readblock(void *p);
};

class AVGCNFile : public AUXFile
{
public:
    AVGCNFile(int t):AUXFile(t){};
private:
    virtual int writeblock(void *p);
    virtual int readblock(void *p);
};

class MDCASKconFile : public AUXFile /* MDCASK configuration file */
{
public:
    MDCASKconFile(int t):AUXFile(t){};
    virtual char * describe();
private:
    virtual int writeblock(void *p);
    virtual int readblock(void *p);
};

class MDCASKinputFile : public AUXFile /* MDCASK input file */
{
public:
    MDCASKinputFile(int t):AUXFile(t){};
    virtual char * describe();
private:
    virtual int writeblock(void *p);
    virtual int readblock(void *p);
};

class PropFile : public AUXFile/* Output property file */
{
public:
    PropFile(int t):AUXFile(t){};
    virtual char * describe();
    virtual int writeentry(void *p);
};


//------------------------------------------------------------------- 
// Define EComplex Class (will be used in Classic Ewald summation)
//-------------------------------------------------------------------
struct EComplex
{
    double Re, Im; //Real and imaginary part (only data member)
public:
    double Norm() { return sqrt(Re*Re+Im*Im); }
    double Norm2() { return Re*Re+Im*Im; }
    
    static void Mul(const EComplex &a, const EComplex &b, EComplex &c)//c=a*b
    {
        double tmp=a.Re*b.Re-a.Im*b.Im;
        c.Im=a.Re*b.Im+a.Im*b.Re; 
        c.Re=tmp;
    }
    static void MulConjugate(const EComplex &a, const EComplex &b, EComplex &c)
    {
        c.Re=a.Re*b.Re+a.Im*b.Im; 
        c.Im=-a.Re*b.Im+a.Im*b.Re;
    }
    EComplex &operator *=(double d)
    {
        Re*=d;
        Im*=d;
        return *this;
    }
    EComplex &operator +=(EComplex &c)
    {
        Re+=c.Re;
        Im+=c.Im;
        return *this;
    }
    EComplex &operator -=(EComplex &c)
    {
        Re-=c.Re;
        Im-=c.Im;
        return *this;
    }
    void ExpI(double theta)
    { 
        Re=cos(theta);
        Im=sin(theta);
    }
    EComplex &operator =(double d) 
    { 
        Re=d;
        Im=0.0;
        return *this;
    }
    friend LOStream &operator <<(LOStream &os, EComplex &c)
    {
        return os << "["<<c.Re<<"(+)"<<c.Im<<"(I)]";
    }
};

/*-------------------------------------------------------------------
 * Define Complex Array Type (will be used in PME)
 *------------------------------------------------------------------- */
typedef double ComplexType[2];

/* MD simulation framework */
class MDFrame : public Organizer
{
public:
    int _SAVEMEMORY;    /* controls how much memory is allocated per atom */
    int _NP;            /* total number of atoms */
    int _NP0, _NIMAGES; /* number of original and image atoms */
    int _NPfree,_NPfixed; /* number of free and fixe atoms */

    class Vector3 *_R, *_SR;   /* real and scaled coordinates of atoms */
    class Vector3 *_SRA; /* running average coordinates of atoms in scaled coordinate */
    class Vector3 *_VR, *_VSR; /* real and scaled velocities of atoms */
    class Vector3 *_F, *_F0;   /* force on atoms */
    class Vector3 *_R0;        /* place to save old coordinates of atoms */
    class Vector3 *_Fext;      /* external force on atoms */
    int _ENABLE_FEXT;

    int *fixed;                /* 1 if atom is fixed */
    int *image;                /* image atom */
    int *species;              /* species type of atom */
    int *group;                /* group ID */
    int nspecies;              /* number of atom species */
    char element[MAXSPECIES][10]; /* species names */

    class Matrix33 _H, _H0;     /* simulation box vectors and copy */
    class Matrix33 _VIRIAL, _VH;/* Virial stress, and velocity of H */
    class Matrix33 *_VIRIAL_IND;/* Virial stress contribution from each atom */
    double _EPOT, *_EPOT_IND;   /* total potential energy and contribution per atom */
    double *_EPOT_RMV;          /* potential energy change if atom is removed */
    double _EPOT_BOX;           /* potential energy of the box */
    double _EPOT0;              /* potential energy from previous step */
    double _EPOT0_ext;          /* reference potential energy from boundary atoms */
    double *_TOPOL;             /* atomic environment, e.g. central symmetry parameter */
    
    /* Molecular Dynamics parameters */
    double _ESTRAIN;             /* box strain energy (in eV) */
    double _OMEGA;               /* box volume (in A^3) */
    double _KATOM, _KBOX;        /* kinetic energies (in eV) */
    double _ATOMMASS[MAXSPECIES], _WALLMASS; /* in (g/mol) */
    double _ATOMCHARGE[MAXSPECIES]; /* in electron */
    double _TIMESTEP;            /* in picosecond (ps) */
    double _TDES,_T;             /* desired and current temperature */
    int _DOUBLET;                /* 1: if initialize T double desired temperature */
    double _ATOMTCPL,_BOXTCPL;   /* temperature coupling constants */
    double _BOXDAMP;             /* damping coefficient for box PR method to prevent oscillation */
    double zeta,zetav,zetaa,vt2; /* nose hoover thermal stat */
    double _HELM,_HELMP;         /* Helmholtz free energy and a conserved var */
//    int usenosehoover;           /* 1: if use nose hoover thermostat */
//    int usescalevelocity;        /* 1: if scale velocity for temperature control */
    int applyshear;              /* apply constant shear deformation */
    double extforce[100];        /* external forces on different atom groups */
    double forcemul;             /* a multiplication factor for extforce */
    Matrix33 _SHEARRATE;         /* shear deformation rate, if applyshear==1 */
    int runavgposition;
    
    /* Neighbor list */
    double _RLIST,_SKIN;        /* potential cut off radius (in A) */
    int _NIC, _current_NIC;     /* maximum allowable and actual number of atoms per cell */
    int _NNM, _current_NNM;     /* maximum allowable and actual number of neighbors per atom */
    int *nn;                    /* nn[i]: number of neighbors of atom i */
    int **nindex;               /* nindex[i][j]: the jth neighbor of atom i */

    int *link_head, *link_list; /* link list */
    
    int ****celllist;
    char *cell_mem, *nindex_mem;
    bool firsttime;

    double rc_plot;             /* another neighbor list for visualization */
    int _NNM_plot;
    int *nn_plot;
    int **nindex_plot;
    char *nindex_plot_mem;

    /* Numerical integrator */
    Vector3 *s2, *s3, *s4, *s5, *s2real; /* Gear 6th predictor-corrector */
    class Matrix33 h2,h3,h4,h5,h2real;
    double zeta2,zeta3,zeta4,zeta5;      /* for Nose-Hoover thermal stat */
    int npold;
    int NHChainLen;                      /* Length Nose-Hoover Chain */
    double zetaNHC[MAXNHCLEN], zetaNHCv[MAXNHCLEN], zetaNHCa[MAXNHCLEN],
           zetaNHC2[MAXNHCLEN], zetaNHC3[MAXNHCLEN], zetaNHC4[MAXNHCLEN], zetaNHC5[MAXNHCLEN];
    double NHMass[MAXNHCLEN];            /* Masses of Nose-Hoover Chain */

                                        /* Nose Hoover chain for box */
//    double zetaBNHC[MAXNHCLEN], zetaBNHCv[MAXNHCLEN], zetaBNHCa[MAXNHCLEN],
//           zetaBNHC2[MAXNHCLEN], zetaBNHC3[MAXNHCLEN], zetaBNHC4[MAXNHCLEN], zetaBNHC5[MAXNHCLEN];    
//    double BNHMass[MAXNHCLEN];         
 
    /* Internal and applied stresses */
    class Matrix33 _PSTRESS;    /* stress due to momenta (in eV/A^3) */
    class Matrix33 _TOTSTRESS;  /* total stress (_PSTRESS+_VIRIAL) */
    class Matrix33 _TOTSTRESSinMPa; /* total stress in MPa */
    class Matrix33 _EXTSTRESS;  /* external applied stress */
    double _TOTPRESSURE;
    double _EXTSTRESSMUL;       /* multiplier of applied stress */
    double _EXTPRESSADD;        /* addition of hydrostatic pressure, not affected by ESTSTRESSMUL */
    double _VACUUMRATIO;        /* stress correction factor when medium is smaller than simulation box */
    class Matrix33 _DVIRIAL_Dexx; /* for computing elastic constants at finite temperature */
    double _Dexx;
    
    /* Conjugate gradient relaxation */   
    int conj_itmax,conj_fevalmax,conj_fixbox,conj_fixshape,
        conj_fixatoms,conj_step,conj_ier,conj_fixdir;
    class Matrix33 conj_fixboxvec;
    double conj_dfpred,conj_f,conj_ftol,conj_g2res;
    double *conj_space;
    double _HPRECOND;                     /* preconditioner for H */
    double _XPRECOND,_YPRECOND,_ZPRECOND; /* preconditioner for x,y,z */
    class Matrix33 _SIGMA;                /* stress term in Parilleno-Raman */
    class Matrix33 _GH;                   /* Energy derivative dE_pot / dH */

    /* Conjugate gradient relaxation with constraints */
    double constrainS, constrainS_inst, constrain_dist, constrainF;
    int constrainatoms[MAXCONSTRAINATOMS];

    /* constrained MD*/
    int _constrainedMD;
    
    /* Molecular Dynamics simulation */
    int totalsteps, equilsteps, curstep, step0, continue_curstep;  /* simulation steps */
    char ensemble_type[100], integrator_type[100]; int implementation_type;
    int algorithm_id;
    
    /* Monte Carlo simulation */
    int MC_atom, MC_accept, MC_accept_tot, MC_switchoffatoms[10];
    int MC_dV_freq, MC_dN_freq;
    double MC_accept_ratio, MC_dVmax, MC_dfrac, MC_P_ext, MC_mu_ext, MC_Lambda3, MC_fracatom;
    Vector3 MC_dr;

    /* Free energy calculations */
    double _LAMBDA, _LAMBDA0, _LAMBDA1, dEdlambda, dlambdadt, _ECOH, _WTOT, _WAVG;
    int _ENABLE_SWITCH, _SWITCHFREQ, _RANDSEED, _CHAINLENGTH, _SWITCHFUNC, _REFPOTENTIAL;
    int heslen, *hesind; double *hes; /* Hessian matrix */
    double _ECSPRING;                 /* spring constant of Einstein crystal */
    double _I12_RC, _I12_SIGMA, _I12_EPSILON; /* Inverse r12 potential */
    double _GAUSS_RC, _GAUSS_SIGMA, _GAUSS_EPSILON; /* Gaussian potential */

    /* Nudged-Elastic-Band (NEB) method for saddle point search */
    class Vector3 *_SR1, *_SR2, *_SR3;
    class Vector3 dscom12;      /* center of mass difference for constrained atoms between _SR1 and _SR1 */
    class Vector3 _COM;         /* center of mass */
    class Vector3 **_Rc, **_Fc; /* atom positions and forces along a chain of states */
    double *_nebinterp;
    int nebspec[10];
    int readparallelnebcn;
    double annealspec[10];

    /* Forward Flux Sampling parameters (FFS) (Seunghwa Ryu) */ 
    double *_QLM, *_QL;         /* Local order parameter for crystaline structure */
    int    *_Lam_array, *_Lam_check; /* lambda arrays for FFS */
    int    L_for_QLM;           /* angular momentum number for QLM */
    int    N_solid_P;           /* total number of solid-like particle */
    int    N_skin_P;            /* total number of skin particle */
    int    N_lgst_cluster;      /* number of particles in largest cluster */
    double Rc_for_QLM;          /* cut-off radius for QLM */
    double QLM_cutoff;          /* order param cutoff for solid-like particle*/
    int    *QLM_mark;           /* DFS search algorithm mark 1 if searched, 0 if not */
    int    *QLM_solid;          /* equals 1 if solid, 0 if not */
    int    N_cluster_temp;      /* temporary memory for cluster size counting */
    Vector3 Principal_Inertia;  /* Inertia for lgst_cluster */
    double N_lgst_index;        /* Special TOPOL indexing for N_lgst_cluster */
    double N_lgst_skin;         /* total number of skin atom for lgst_cluster*/
    double N_lgst_skin_index;   /* Special TOPOL indexing for lgst skin */
  
    int    wSKIN;               /* 1 for skin inclusion for cluster, 0 for not */  
    int    saveFFScn;           /* configuration save if saveFFScn=1 */
    int    FFSoption;           /* 0 for q/T, 1 for P(i+1/i) */
    int    FFSpruning;          /* 0 for time-limited-run, 1 for pruning */
    double    FFScn_weight;     /* default = 1, pruning factor */
    double    FFS_Pp;           /* Pruning Probability */
    int    FFS0_check;          /* check key for FFS0 run */
    int    FFS_i;               /* FFS i th order parameter */
    int    FFS_i_plus;          /* FFS (i+1)th order parameter */
    int    FFS_i_minus;         /* FFS (i-1)th order parameter */
    
    /* Ewald Summation for Coulomb interaction */
    Vector3 *_F_Ewald, *_F_Ewald_Real, *_F_Ewald_Rec;
    double *_EPOT_IND_Ewald, *_EPOT_IND_Ewald_Real, *_EPOT_IND_Ewald_Rec;
    Matrix33 _VIRIAL_Ewald, _VIRIAL_Ewald_Real, _VIRIAL_Ewald_Rec;
    int Ewald_CE_or_PME;    /*0: CE (classical Ewald), 1: PME (particle mesh Ewald) */
    int Ewald_option_Alpha; /*different options to determine Alpha (spread of Gaussian charge) */
    double Ewald_Rcut;      /* real space cut-off distance for Ewald summation */
    double Ewald_precision; /* cut-off error in both real and reciprocal space */
    double Ewald_time_ratio;/* computational time ratio between real and reciprocal space */
    double _EPOT_Ewald, _EPOT_Ewald_Real, _EPOT_Ewald_Rec, _EPOT_Ewald_Self;
    double Ewald_Alpha, Ewald_qSQRsum;
    double Ewald_H_L1, Ewald_H_L2, Ewald_H_L3, Ewald_cell_V; /* simulation box size */
    
    int CE_MKV, CE_nKV;         /* max and actual number of K-vectors in list */
    int CE_Mka, CE_Mkb, CE_Mkc; /* Largest repeat reciprocal distance in each direction */
    double CE_Kc;               /* cut-off radius in reciprocal space */
    
    struct K_rec { double kx, ky, kz, k2; int a,b,c; } *CE_KV; /* Actual K-vectors */
    EComplex *CE_sc_axby, *CE_sckk;  /* structural factors exp(iK(rx+ry)), qj*exp(iK(rx+ry+rz)) */
    EComplex **CE_scx, **CE_scy, **CE_scz; /* sin-cos table for computing structural factor */
    
    clock_t CE_sREAL, CE_sREC, CE_eREAL, CE_eREC; /* time information for CE calculation */
    double  CE_iREAL, CE_iREC, CE_tREAL, CE_tREC;
    
    /* Particle Mesh Ewald Summation (Hark Lee, 2007/10/18) */
    int PME_bsp_n;                       /* order of B-spline interpolation */

    int PME_K1, PME_K2, PME_K3, PME_K3S; /* number of mesh points */
    int PME_K23, PME_K123, PME_K23S, PME_K123S;
    double PME_K1d, PME_K2d, PME_K3d, PME_K123d;
    
    Vector3 *_UR;                        /* Scaled atom positions */

    double *PME_B;                       /* B array */
    double *PME_C;                       /* C array */
    double *PME_BC;                      /* B*C array */
    double *PME_Q;                       /* Q array */
    double *PME_CONV;                    /* convolution of BC and IQ */
    ComplexType *PME_IQ;                 /* inverse Q array */

    ComplexType PME_R_bsp_Bm;            /* entries of B */
    double *PME_bsp_Mk;                  /* to calculate B entries */
    double *PME_bsp_fac;                 /* factorial array */

    int *PME_m1s, *PME_m2s, *PME_m3s;    /* k-space start indices */
    int *PME_m1, *PME_m2, *PME_m3;       /* k-space indices */
    int *PME_m1mod, *PME_m2mod, *PME_m3mod; /* actual indices after PBC */

    double *PME_x1, *PME_x2, *PME_x3;    /* (u-k) after PBC */
    double *PME_MU1, *PME_MU2, *PME_MU3; /* Mn(u) */
    double *PME_d1MU1, *PME_d1MU2, *PME_d1MU3; /* Mn-1(u) */
    double *PME_d2MU1, *PME_d2MU2, *PME_d2MU3; /* Mn-1(u-1) */

    fftw_plan PME_fft_plan1, PME_fft_plan2; /* fftw plans */
    
    clock_t PME_sREC, PME_eREC;          /* time information for PME calculation */
    double  PME_iREC, PME_tREC;
    
    /* Configuration manipulation */
    class Vector3 *storedr;         /* atomic displacements */
    class Matrix33 dH;              /* change of matrix H (cell size) */
    char  crystalstructure[30];     /* crystal structure for makecrystal() */
    double latticeconst[3];         /* lattice constant  for makecrystal() */
    double latticesize[3][4];       /* lattice size for makecrystal() */
    int _TORSIONSIM;
    double torquespec[3];           /* torsion specification */
    double _TORQUE, *_TORQUE_IND;
    int _BENDSIM;
    double bendspec[4];             /* bending specification */
    double _BENDMOMENT, *_BENDMOMENT_IND;
    Vector3 _P_COM, _F_COM;         /* center of mass momenta */
    double _PTHETA_COM; 

    /* File input and output */
    double input[20000];                    /* general input variable in script file */
    double output_dat[2000];                /* general output data array */
    char   output_str[10000];               /* general output string */
    char   output_fmt[10000];               /* output format */
    class CNFile initcn,intercn,FFScn,finalcn;/* configuration file objects */
    class AVGCNFile avgcn;/* file object that contains running averaged position */
    class PropFile pf;                      /* property file objects */
    char incnfile[200], finalcnfile[200];   /* configuration file name */
    char intercnfile[200];                  /* configuration file name */
    char FFScnfile[200];
    char intercfgfile[200];                 /* atomeye cfg file name */
    char outpropfile[200];                  /* property file name */
    char myname[200];                       /* name of the simulation */
    char potfile[100];                      /* name of potential file */    
    char command[MAXCMDLEN]; int ncom;           /* shell command */
    int  savecn, savecnfreq;                /* frequency of saving cn files */
    int  savecfg, savecfgfreq, savecfgcount;/* frequency of saving cfg files */
    int  saveprop, savepropfreq;            /* frequency of saving prop files */
    int  printfreq;                         /* frequency of printing on scren */
    int  filecounter;                       /* number of current file */
    int  FFSfilecounter;                    /* number of current FFS file */
    int allocmultiple;                      /* allocate more memory than current number of atoms */
    int writevelocity;                      /* write velocity into cn file */
    int writeall;                           /* write "sx,sy,sz,vx,vy,vz,pot,fixed" into cn file */
    int zipfiles;                           /* zip output files */
    int fixedatomenergypartition;           /* 1: disregard energy contribution from fixed atoms */    

    /* Interface to Fortran code (Relax) */
    char fortranpath[100], fortranexe[100];

    /* Interface to AtomEye and other viewers */
    char atomeyepath[100], atomeyeexe[100];
    int atomeyerepeat[4];

    /* Interface to VASP */
    class POSCARFile initposcar,finalposcar;

    /* Interface to LAMMPS */
    class LAMMPSFile initlammps,finallammps;
    
    /* Interface to MDCASK */
    class MDCASKconFile initmdcaskcon,finalmdcaskcon;
    class MDCASKinputFile initmdcaskinput,finalmdcaskinput;
    char MDCASKpot[50];
        
    /* Visualization */
    YWindow *win;                   /* Display window object */

    int win_width, win_height;      /* size of the window */
    double rotateangles[5];         /* rotation angles and scaling parameter */
    int plotfreq;                   /* frequency of re-plotting */
    double atomradius[MAXSPECIES];  /* size of atoms */
    double bondradius, bondlength;  /* thickness and length of bonds */

    char atomcolor[MAXSPECIES][30]; /* color of atom species */
    char bondcolor[30];             /* color of bonds */    
    char fixatomcolor[30];          /* color of fixed atoms, fixed[i]=1 */
    char highlightcolor[30];        /* color of highlighted atoms */
    char backgroundcolor[30];       /* background color of window */
    char colornames[MAXCOLORS][30]; /* names of allocated colors */
    unsigned colors[MAXCOLORS+15];  /* value of allocated colors */

    int plot_highlight_atoms[10000];/* indices of atoms to be highlighed */    
    double plot_limits[10];         /* set x,y,z, limits of plotting regime */
    int    plot_atom_info;          /* print information when atom is clicked */
    int    plot_map_pbc;            /* map all atoms in primary cell when plotting */
    double plot_color_windows[100]; /* only plot atoms whose properties fall into specified windows */
    double plot_color_bar[10];      /* how to color atoms according to their properties */
    int    plot_color_axis;          /* which atomic property specifies color */

    int autowritegiffreq;           /* frequency of outputing .GIF graphics file */
    double *color_ind;              /* property linked to atom color */
    int NCS;                        /* number of neighbors for central symmetry parameter comp. */
    int coloratoms;                 /* determine how to color atoms 
                                   0 = color by species, 1 = color by groups */

    /* Constructor (set initial values of member variables) */
    MDFrame(): _SAVEMEMORY(0),_NP(0),_NP0(0),_NIMAGES(0),_NPfree(0),_NPfixed(0),

               _R(0),_SR(0),_SRA(0),_VR(0),_VSR(0),_F(0),_R0(0),_Fext(0),_ENABLE_FEXT(0),

               fixed(0),species(0),group(0),nspecies(1),
               _EPOT(0),_EPOT_IND(0),_EPOT_RMV(0),_EPOT_BOX(0),_EPOT0(0),_EPOT0_ext(0),_TOPOL(0),

               /* Molecular Dynamics parameters */
               _ESTRAIN(0),_KATOM(0),_KBOX(0),_WALLMASS(1),
               _TIMESTEP(1e-3),_TDES(0),_T(0),_DOUBLET(0),
               _ATOMTCPL(100),_BOXTCPL(100),_BOXDAMP(0),
               zeta(0),zetav(0),zetaa(0),vt2(0.25),_HELM(0),_HELMP(0),
               //usenosehoover(0),usescalevelocity(0),
               applyshear(0),
               forcemul(1), runavgposition(0),

               /* Neighbor list */
               _RLIST(0),_SKIN(0),_NIC(100),_current_NIC(0),_NNM(60),_current_NNM(0),nn(0),nindex(0),
               link_head(0),link_list(0),
               celllist(0), cell_mem(0), nindex_mem(0), firsttime(true),
               rc_plot(0), _NNM_plot(4), nn_plot(0), nindex_plot(0), nindex_plot_mem(0),
               
               /* Numerical integrator */
               s2(0),s3(0),s4(0),s5(0),s2real(0), /* Gear 6th predictor-corrector */
               zeta2(0),zeta3(0),zeta4(0),zeta5(0),
               npold(0),

               /* Nose Hoover chain */
               NHChainLen(2),

               /* Internal and applied stresses */
               _TOTPRESSURE(0),_EXTSTRESSMUL(1),_EXTPRESSADD(0),_VACUUMRATIO(0),_Dexx(0),

               /* Conjugate gradient relaxation */
               conj_itmax(1000),conj_fevalmax(10000),conj_fixbox(1),
               conj_fixshape(0),conj_fixatoms(0),conj_step(0),conj_ier(0),conj_fixdir(0),
               conj_dfpred(1e-3),conj_f(0),conj_ftol(1e-8),conj_g2res(0),conj_space(0),
               _HPRECOND(1),_XPRECOND(1),_YPRECOND(1),_ZPRECOND(1),

               /* Conjugate gradient relxation with constrains */
               constrainS(0),constrain_dist(0),
               
               /* constrained MD */
               _constrainedMD(0),
               
               /* Molecular Dynamics simulation steps */
               totalsteps(100), equilsteps(0), curstep(0), step0(0), continue_curstep(0),
               implementation_type(0), algorithm_id(0),

               /* Monte Carlo simulation */
               MC_atom(0),MC_accept(0),MC_accept_tot(0), 
               MC_dV_freq(1), MC_dN_freq(1), MC_accept_ratio(0), MC_dVmax(0.1), MC_dfrac(0.1),
               MC_P_ext(0), MC_mu_ext(0), MC_Lambda3(0), MC_fracatom(1),

               /* Free energy calculations */
               _LAMBDA1(1),dEdlambda(0),dlambdadt(0),_ECOH(0),_WTOT(0),_WAVG(0),
               _ENABLE_SWITCH(0), _SWITCHFREQ(1),
               _RANDSEED(12345),_CHAINLENGTH(0),_SWITCHFUNC(0),_REFPOTENTIAL(0),
               heslen(0), hesind(0), hes(0),

               /* Switching potentials */
               _ECSPRING(0),
               _I12_RC(0), _I12_SIGMA(0), _I12_EPSILON(0),
               _GAUSS_RC(0), _GAUSS_SIGMA(0), _GAUSS_EPSILON(0),
               
               /* Nudged-Elastic-Band (NEB) method */
               _SR1(0),_SR2(0),_SR3(0),_Rc(0),_Fc(0),_nebinterp(0), //_Fnebinterp(0),
               readparallelnebcn(0),

               /* Forward Flux Sampling (FFS) parameters (Seunghwa Ryu) */ 
               _QLM(0),_QL(0),_Lam_array(0),_Lam_check(0), L_for_QLM(0),
               N_solid_P(0),N_skin_P(0),N_lgst_cluster(0),Rc_for_QLM(0),
               QLM_cutoff(0),QLM_mark(0), QLM_solid(0), N_cluster_temp(0), 
               N_lgst_index(0), N_lgst_skin(0), N_lgst_skin_index(0), wSKIN(0),
               saveFFScn(0), FFSoption(0), FFSpruning(0), FFScn_weight(1), 
               FFS_Pp(0), FFS0_check(1), FFS_i(0), FFS_i_plus(0), 
               FFS_i_minus(0), 
               
               /* Ewald Summation for Coulomb interaction */
               _F_Ewald(0), _F_Ewald_Real(0), _F_Ewald_Rec(0),
               _EPOT_IND_Ewald(0), _EPOT_IND_Ewald_Real(0), _EPOT_IND_Ewald_Rec(0),
               Ewald_CE_or_PME(0),
               Ewald_option_Alpha(0),
               Ewald_Rcut(9),
               Ewald_precision(4.2),
               Ewald_time_ratio(3.6),
               _EPOT_Ewald(0), _EPOT_Ewald_Real(0),
               _EPOT_Ewald_Rec(0), _EPOT_Ewald_Self(0),               
               Ewald_Alpha(0), Ewald_qSQRsum(0),
               Ewald_H_L1(0), Ewald_H_L2(0), Ewald_H_L3(0), Ewald_cell_V(0), 
               
               CE_MKV(0), CE_nKV(0),
               CE_Mka(0), CE_Mkb(0), CE_Mkc(0),
               CE_Kc(0),
               CE_KV(0),
               CE_sc_axby(0), CE_sckk(0),
               CE_scx(0), CE_scy(0), CE_scz(0),
               CE_sREAL(0), CE_sREC(0), CE_eREAL(0), CE_eREC(0),
               CE_iREAL(0), CE_iREC(0), CE_tREAL(0), CE_tREC(0),    
               
               /* Particle Mesh Ewald Summation (Hark Lee, 2007/10/18) */ 
               PME_bsp_n(8),
               PME_K1(64), PME_K2(64), PME_K3(64), PME_K3S(0),
               PME_K23(0), PME_K123(0), PME_K23S(0), PME_K123S(0),
               PME_K1d(0), PME_K2d(0), PME_K3d(0), PME_K123d(0),
               _UR(0),

               PME_B(0), PME_C(0), PME_BC(0), PME_Q(0), PME_CONV(0), PME_IQ(0),
               PME_bsp_Mk(0), PME_bsp_fac(0),
               PME_m1s(0), PME_m2s(0), PME_m3s(0),
               PME_m1(0), PME_m2(0), PME_m3(0),
               PME_m1mod(0), PME_m2mod(0), PME_m3mod(0),

               PME_x1(0), PME_x2(0), PME_x3(0),
               PME_MU1(0), PME_MU2(0), PME_MU3(0),
               PME_d1MU1(0), PME_d1MU2(0), PME_d1MU3(0),
               PME_d2MU1(0), PME_d2MU2(0), PME_d2MU3(0),

               PME_sREC(0), PME_eREC(0), PME_iREC(0), PME_tREC(0),
               
               /* Configuration manipulation */
               storedr(0),

               /* Torsion simulation */
               _TORSIONSIM(0),_TORQUE(0),

               /* Bending simulation */
               _BENDSIM(0), _BENDMOMENT(0),
        
               /* Center of mass momenta */
               _PTHETA_COM(0),               

               /* File input and output */
               initcn(AUXFile::BLOCK),intercn(AUXFile::SERIES),
               FFScn(AUXFile::SERIES),
               finalcn(AUXFile::BLOCK),avgcn(AUXFile::BLOCK),
               pf(AUXFile::ENTRIES),
               ncom(0),savecn(0),savecnfreq(100),savecfg(0),
               savecfgfreq(100),savecfgcount(1),
               saveprop(0),savepropfreq(100),
               printfreq(100),filecounter(1),FFSfilecounter(1),
               allocmultiple(1),writevelocity(0),writeall(0),
               zipfiles(0),fixedatomenergypartition(0),

               /* Interface to VASP */
               initposcar(AUXFile::BLOCK),finalposcar(AUXFile::BLOCK),

               /* Interface to LAMMPS */
               initlammps(AUXFile::BLOCK),finallammps(AUXFile::BLOCK),

               /* Interface to MDCASK */
               initmdcaskcon(AUXFile::BLOCK),finalmdcaskcon(AUXFile::BLOCK),
               initmdcaskinput(AUXFile::BLOCK),finalmdcaskinput(AUXFile::BLOCK),

               /* Visualization */
               win(0),win_width(350),win_height(350),
               plotfreq(1),bondradius(0.1),bondlength(0),plot_atom_info(1),plot_map_pbc(0),
               plot_color_axis(0),autowritegiffreq(0),color_ind(0),NCS(12),coloratoms(0)
               
    {
        /* Input/output control */
        input[0]=0;
        output_str[0]=0;
        output_fmt[0]=0;
        _ATOMMASS[0]=1;
        
        /* Configuration manipulation */
        latticeconst[0]=latticeconst[1]=latticeconst[2]=1.0;
        sprintf(crystalstructure,"simple-cubic");

        /* Workspace for Monte Carlo simulation */
        MC_switchoffatoms[0]=0;

        /* Plot settings */
        plot_highlight_atoms[0]=0;
        plot_limits[0]=0;
        plot_color_windows[0]=0;
        plot_color_bar[0]=0;

        
    };

    virtual ~MDFrame() { delete(win); }
            
    /* Parser */    
    virtual void initparser();
    virtual int exec(char *name);            
    virtual void initvars();
    void runcommand();        /* run shell command */
#ifdef _USETCL    
    virtual int Tcl_AppInit(Tcl_Interp *interp);
#endif
    
    /* Coordinate transformation */
    void SHtoR();
    void RHtoS();
    void RtoR0();
    void R0toR();
    void clearR0(); /* clear R0 so that neighbor list will be reconstructed */
    bool Bond(int I, int J) const;
    
    virtual void Alloc();

    /* Potentials */    
    virtual void call_potential(); /* by default calls potential but in mdparallel calls potential_parallel */
    virtual void potential();
    virtual void potential_energyonly();
    virtual double potential_energyonly(int iatom);
    void ECpotential();            /* Harmonic potential of Einstein crystal */
    void ECpotential_energyonly();
    double ECpotential_energyonly(int iatom);
    void HApotential();
    void HApotential_energyonly(); /* Harmonic potential based on Hessian matrix */
    void I12potential_energyonly(); /* Inverse r-12 potential as liquid reference */
    void I12potential();
    void GAUSSpotential();          /* Gaussian potential as liquid reference */
    void SWITCHpotential_energyonly(double lambda); /* for adiabatic switching */
    void SWITCHpotential(double lambda); 
    virtual void SWITCHpotential_user(double lambda); /* user defined function for switching */
    virtual void SWITCHpotential_user_energyonly(double lambda);
    double SWITCH_Find_Lambda();
    void eval_insertparticle();      /* compute potential energy of test particle */
    void eval_removeparticle();      /* compute potential energy of test particle */

    /* Neighbor list */
    virtual void NbrList_reconstruct(int iatom=-1);  
    void NbrList_reconstruct_use_cell_list(int iatom=-1); /* use cell list with Verlist */
    void NbrList_reconstruct_use_link_list(int iatom=-1); /* use link list with Verlist */
    
    void NbrList_print();        /* for debug use */
    void NbrList_refresh();
    bool NbrList_needrefresh();
    void NbrList_init(int, int);
    void NbrList_free();

    void NbrList_initlinklist(int);
    void NbrList_freelinklist();
    
    void NbrList_initcell(int,int,int,int);
    void NbrList_freecell();

    void NbrListPlot_reconstruct();
    void NbrListPlot_init(int,int);
#define refreshneighborlist NbrList_refresh
    
    /* Numerical integrator */
    void NVE_Gear6();
    void NVT_Gear6();
    void NVTC_Gear6();
    void NPH_Gear6();
    void NPT_Gear6();
    void NPTC_Gear6();
    void Gear6_init();
    void VVerlet_Get_s2();
    void VVerlet_Get_h2();
    void VVerlet_Get_Tinst();
    void VVerlet_NHCINT();

    void NVE_VVerlet(int);
    void NVT_VVerlet_Implicit(int);
    void NVT_VVerlet_Explicit_1(int);
    void NVT_VVerlet_Explicit_2(int);
    void NVTC_VVerlet_Explicit(int);
    void NPH_VVerlet(int);
    void NPH_VVerlet_Explicit_1(int);
    /* void NPH_VVerlet_Explicit_2(int); */ /* not working yet */
    void NPT_VVerlet_Implicit(int);
    void NPT_VVerlet_Explicit_1(int);
    void NPT_VVerlet_Explicit_2(int);

    /* Conjugate gradient relaxation */
    int applyconstraint();
    int relax();           /* CG relax */
    int constrainedrelax();/* CG relax with constraints */
    void relax_by_group(double *);
    static void potential_wrapper(int *n,double x[],double *f,double df[]);
    static void potential_wrapper_c(int n,double x[],double *f,double df[]);
    static void cstr_potential_wrapper(int *n,double x[],double *f,double df[]);
    static void cstr_potential_wrapper_c(int n,double x[],double *f,double df[]);
    void strelax();         /* steepest descent relaxation */
    
    /* Molecular Dynamics simulation */
    virtual void integrator();       /* MD integrator */

    virtual void run();              /* run the simulation */
    virtual void calcprop();         /* calculate physical properties */
    virtual void calcoutput();       /* calculate output string */
    void printResult();
    void calDVIRIALDexx();           /* calculate strain derivative of stress */
    void calcentralsymmetry();       /* calculate central symmetry deviation */
    void caldisregistry();           /* for slip (dislocation) visualization */
    void calHessian();               /* numerical calculation of Hessian matrix */
    void readHessian();              /* read Hessian matrix from file */
    void calModeHessian();           /* calculate Hessian matrix for a set of deformation modes */
    void calphonondisp();            /* calculate phonon dispersion relation */
    void calmisfit();                /* calculate misfit potential */
    void calmisfit2();               /* calculate misfit potential */
    void testpotential();            /* test consistency between force and energy */
    void findcore();                 /* find center of mass of defect core atoms */
    void initvelocity();
    void randomposition();
    //void velscaler(); /* disable velocity scaler */
    void perturbevelocity();         /* gaussian perturbation */
    void MCperturbevelocity();       /* shooting for MC TPS for NVE, NVT */
    void multiplyvelocity();
    void step();
    void eval();
    void multieval();                /* for testing/timing purposes */

    /* Monte Carlo simulation */
    int runMC();                    /* Metropolis Monte Carlo */
    int runTAMC();                  /* temperature accelerated Monte Carlo */
    int MCstep();
    int MCstep_moveatom();
    int MCstep_dV();
    int MCstep_dN();
    int MCstep_fracatom();
    int MCstep_moveatom_fracatom();
    int MCstep_dV_fracatom();
    int MCstep_dN_fracatom();    
    int runMCSWITCH();
    int MCSWITCHstep(double lambda);
    int runMDSWITCH();
    
    /* Nudged-Elastic-Band relaxation */
    void AllocChain();
    void caldscom12();
    void initRchain();
    int readRchain();
    int writeRchain();
#ifdef _STRINGMETHOD
    void runstringmethod();
#endif
    void nebrelax();   
    void annealpath();       /* simulated annealing of paths */
    void cutpath();
    int statedistance();     /* compute distance between two states _SR1 and _SR */
    int interpCN(double s); /* 2007/3/6 Wei Cai */
    void copyRchaintoCN(int j);
    void copyCNtoRchain(int j);
    void moveRchain(int, int);

    /* Forward Flux Sampling (FFS) functions */
    void calcrystalorder();
    void allocQLM();
    void assign_Lam();
    void DFS(int,int,double);
    void DFS1(int,int,double,double); 
    int FFSprocess();

    /* Ewald Summation for Coulomb interaction */
    void Ewald_init();
    void CE();
    void CE_Real();
    void CE_Rec();
    void CE_init();
    void CE_clear();
    void CE_enumKV();
    void CE_fillSinCosTables();
    
    /* Particle Mesh Ewald Summation (Hark Lee, 2007/10/18) */
    void PME();
    void PME_Rec();
    void PME_init();
    void PME_clear();
    void PME_cal_bsp_fac();
    void PME_cal_bsp_Mk();
    double PME_cal_bsp_Mu(int n, double u);
    void PME_cal_bsp_Bm(int mi, int Ki);
    double PME_cal_bsp_Cm(int m1, int m2, int m3);
    
    /* Configuration manipulation */
    void makecrystal();         /* create perfect crystal structure */
    void makecut();             /* create cut surfaces in lattice */
    void makedipole();          /* create dislocation dipole structure */
    void makedislocation();     /* create an arbitrarily oriented dislocation line */
    void makecylinder();        /* cut a cylinder out of the crystal */
    void setfixbufferatoms();   /* fix atoms around a cylinder */
    void makedislcylinder();    /* create a single dislocation in a cylinder */
    
    void makedisloop();         /* create dislocation loop */
    void calloopdisp            /* compute the displacement field */
         (Vector3 *,double,double,Vector3 *,Vector3 *,Vector3 *);
    void makegrainboundary();   /* create grain boundary */
    void makewave();            /* initiate wave form */
    void maketorquehandle();    /* prepare for torque application */
    void removetorquehandle();  /* restore usual PBC */
    void makebendhandle();      /* prepare for bend application */
    void copytorqueatoms();     /* copy end atoms for torsion simulation */
    void copybendatoms();       /* copy end atoms for bending simulation */
    
    void scaleH();              /* hydrostatic scale of simulation box */
    void setH();
    void saveH();
    void restoreH();            /* restore H from H0 */
    void scaleVel();            /* scale velocity */
    
    void shiftbox();            /* shift PBC box */
    void redefinepbc();
    void applystrain();         /* apply homogeneous strain */
    void extendbox();           /* make multiple copies of the box*/

    void cutslice();            /* cut a portion from box, inverse of extendbox */
    void splicecn();            /* put two configurations together */
    void cutpastecn();          /* cut and paste a region of the config */

    void setconfig1();          /* copy current _SR to _SR1 */
    void setconfig2();          /* copy current _SR to _SR2 */
    void switchconfig();        /* R - R1 + R2 */
    void copytoconfig1(int);
    void copytoconfig2(int);
    void replacefreeatom();     /* replace free (not fixed) atoms from another config */
    void relabelatom();
    
    void moveatom();
    void movegroup();
    void setgroupcomvel();      /* set the center of mass velocity of specified groups */
    void addtorque();
    void addbend();
    void printatoms_in_sphere();/* print the indices of all atoms in a sphere */
    void pbcshiftatom();        /* shift atoms in PBC box */
    
    void fixatoms_by_ID();      /* fix a set of atoms whose ID are specified */
    void fixatoms_by_position();/* fix all atoms whose position falls within a specified regime */
    void fixatoms_by_group(double *);
    void fixatoms_by_species(double *);
    void fixatoms_by_pos_topol();/* fix atoms whose position and TOPOL value are in a specified regime */
    void fixallatoms();
    void freeallatoms();
    void reversefixedatoms();
    void constrain_fixedatoms();
    void fix_constrainedatoms();

    void setfixedatomsspecies();
    void setfixedatomsgroup();
    void reversespecies();      /* species: 0->1 , 1->0 */
    void movefixedatoms();      /* move atoms that are marked as fixed */
    void removefixedatoms();    /* remove atoms that are marked as fixed */
    void markremovefixedatoms();/* mark fixed atoms for removal (fix=-1) */
    void removeellipsoid();     /* remove all atoms within an ellipsoid */    
    void removerectbox();       /* remove all atoms within a parallelpepid */
    void removeoverlapatoms();  /* remove atoms that are too close to each other*/
    void findcenterofmass(class Vector3 *);
    void translate_com();       /* apply rigid-body translation so that COM matches that of config1 */
    void rotate_com();          /* apply rigid-body rotation so that COM matches that of config1 */

    void clearFext();           /* reset Fext array */
    void addFext_to_group();    /* add Fext to atoms belonging to specific group */

    /* File input and output */
    virtual void saveintercn(int);   /* save intermediate cn files */
    int readcn();
    int readPOSCAR();
    int readMDCASK();
    int readMDCASKJAIME();
    int readXYZ();
    int readLAMMPS();
    int convertXDATCAR();
    int setfilecounter();
    int setFFSfilecounter();
    int openintercnfile();
    int openFFScnfile(); 
    int openpropfile();
    int closepropfile();
        
    int writefinalcnfile(int zip=1,bool bg=true);
    int writeavgcnfile(int zip=1,bool bg=true);
    int writePOSCAR();
    int writeLAMMPS();
    int writeMDCASK();
    void writePINYMD(char *fname);
    int writeintercnfile(int zip=1,bool bg=false);
    int writeFFScnfile(int zip=1,bool bg=false);

    /* Interface with Fortran code (Relax) */
    void writefortraninifile(char *fname);
    void fortranrelax();

    /* Interface with AtomEye and other viewers */
    virtual void writeatomeyecfgfile(char *fname);
    void convertCNtoCFG();
    void atomeye();
    void writepovray(char *fname);
    void writeMDCASKXYZ(char *fname);
    void writeRASMOLXYZ(char *fname);
    void readRASMOLXYZ(char *fname);
    void writeatomtv(char *fname);
    
    /* Print out energies of current configuration */
    void writeENERGY(char *fname);
    void writeFORCE(char *fname);
    void writePOSITION(char *fname);
    void GnuPlotHistogram();

    /* Visualization */
    virtual void winplot();        
    virtual void winplot(int);
    virtual void openwin();
    virtual void plot();
    int     openwindow(int w,int h,const char *n);
    void    closewindow();
    void    wintogglepause();
    void    alloccolors();
    void    alloccolorsX();
    void    rotate();
    void    saverot();
    
};


/* TCL driver program */
#ifdef _USETCL

int MD_Cmd(ClientData clientData,Tcl_Interp *interp,int argc, const char *argv[]);
int MD_SetVar(ClientData clientData,Tcl_Interp *interp,int argc, const char *argv[]);
int MD_GetVar(ClientData clientData,Tcl_Interp *interp,int argc, const char *argv[]);
int Tcl_Main_Wrapper(int argc, char *argv[]);

#define Tcl_Parser_Init(sim)                                         \
int Tcl_AppInit(Tcl_Interp *interp)                                  \
{                                                                    \
    return sim.Tcl_AppInit(interp);                                  \
}                                                                    \
                                                                     \
int MD_Cmd(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[])\
{                                                                    \
    FILE *istream, *ostream; int mypipe[2];                          \
                                                                     \
    if (pipe (mypipe))                                               \
    {                                                                \
        fprintf (stderr, "Pipe failed.\n");                          \
        return EXIT_FAILURE;                                         \
    }                                                                \
    istream = fdopen(mypipe[0], "r");                                \
    ostream = fdopen(mypipe[1], "w");                                \
    for(int i=1;i<argc;i++)                                          \
        fprintf(ostream,"%s ",(char *)argv[i]);                      \
    fprintf(ostream,"\n");                                           \
    fclose(ostream);                                                 \
    sim.parse_line(istream);                                         \
    fclose(istream);                                                 \
    return TCL_OK;                                                   \
}                                                                    \
                                                                     \
int MD_SetVar(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[])\
{                                                                    \
    int i, curn;                                                     \
    curn = sim.identify(argv[1]);                                    \
    if(curn<0)                                                       \
    {                                                                \
        WARNING("Unrecognized variable "<<argv[1]);                  \
        return TCL_OK;                                               \
    }                                                                \
    for(i=2;i<argc;i++)                                              \
    {                                                                \
        switch(sim.vartype[curn])                                    \
        {                                                            \
        case(INT): sscanf(argv[i],"%d",(int *)sim.varptr[curn]+i-2+sim.shift); break; \
        case(LONG): sscanf(argv[i],"%ld",(long *)sim.varptr[curn]+i-2+sim.shift); break;\
        case(DOUBLE): sscanf(argv[i],"%lf",(double *)sim.varptr[curn]+i-2+sim.shift); break;\
        case(STRING): strcpy((char *)sim.varptr[curn]+i-2+sim.shift,argv[i]); break;   \
        default: FATAL("unknown vartype ("<<sim.vartype[curn]<<")"); \
        }                                                            \
    }                                                                \
    return TCL_OK;                                                   \
}                                                                    \
                                                                     \
int MD_GetVar(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[])\
{                                                                    \
    int i, curn, offset;  char buf[1000];                            \
    curn = sim.identify(argv[1]);                                    \
    if(curn<0)                                                       \
    {                                                                \
        WARNING("Unrecognized variable "<<argv[1]);                  \
        return TCL_OK;                                               \
    }                                                                \
    for(i=1;i<argc;i++)                                              \
    {                                                                \
        if(i==1) offset=0;                                           \
        else sscanf(argv[i],"%d",&offset);                           \
        if(i>2) sprintf(buf," ");                                    \
        switch(sim.vartype[curn])                                    \
        {                                                            \
        case(INT): sprintf(buf,"%d",*((int *)sim.varptr[curn]+offset+sim.shift)); break; \
        case(LONG): sprintf(buf,"%ld",*((long *)sim.varptr[curn]+offset+sim.shift)); break;\
        case(DOUBLE): sprintf(buf,"%.16g",*((double *)sim.varptr[curn]+offset+sim.shift)); break;\
        case(STRING): sprintf(buf,"%s",(char *)sim.varptr[curn]+offset+sim.shift); break;\
        default: FATAL("unknown vartype ("<<sim.vartype[curn]<<")"); \
        }                                                            \
    }                                                                \
    Tcl_SetResult(interp, buf, TCL_VOLATILE);                        \
    return TCL_OK;                                                   \
}                                                                    

#endif

#endif //_MD_H

---