step_defn.h

/*
 * Copyright (C) 2007 by David J. Hardy.  All rights reserved.
 */

#ifndef STEPDEFN_H
#define STEPDEFN_H

#include "step/step.h"
#include "random/random.h"

#ifdef __cplusplus
extern "C" {
#endif

  /* Boolean values */
  enum { STEP_FALSE = 0, STEP_TRUE = 1 };

  /* indices for kinetic virial tensor */
  enum {
    VIRIAL_XX = 0,
    VIRIAL_XY,
    VIRIAL_XZ,
    VIRIAL_YX,
    VIRIAL_YY,
    VIRIAL_YZ,
    VIRIAL_ZX,
    VIRIAL_ZY,
    VIRIAL_ZZ,
    VIRIAL_LEN
  };

  /* indices for upper triangle of kinetic virial */
  enum {
    VIRIAL_UPPER_XX = 0,
    VIRIAL_UPPER_XY,
    VIRIAL_UPPER_XZ,
    VIRIAL_UPPER_YY,
    VIRIAL_UPPER_YZ,
    VIRIAL_UPPER_ZZ,
    VIRIAL_UPPER_LEN
  };

  /* indices for Step::sham array */
  enum {
    SHADOW2=0, SHADOW4, SHADOW6, SHADOW8,
    SHADOW10, SHADOW12, SHADOW14, SHADOW16,
    SHADOW18, SHADOW20, SHADOW22,  SHADOW24,

    ORDERK, ORDERK1
  };

  /* Drude helper class */
  typedef struct DrudeBond_t {
    double massRatio;
    int32 atomIndex;
    int32 drudeIndex;
  } DrudeBond;


  struct Step_t {
    StepParam *param;       /* shallow copy of user StepParam object */
    StepSystem *system;     /* allocated */
    int32 cleanup_vel;      /* did we alloc vel array? */
    int32 cleanup_force;    /* did we alloc force array? */
    int32 doFixedAtoms;
    int32 doRigidWaters;

    Random random;          /* random number generator */

#if 0
    int32 is_force_valid;
#endif

    double tempkonst;       /* constant needed to compute temperature */
    double *scal_inv_mass;  /* compute array of scaled inverse mass,
                               for efficiency */
    MD_Dvec *half_vel;      /* half-step velocities for kinetic reductions */

    /* buffer for error messages */
    char errmsg[512];

    /* retain method compute and done routines */
    int (*compute)(struct Step_t *, int32 numsteps);
    void (*done)(struct Step_t *);

    /* 
     * Variables for calculating shadow Hamiltonian using
     * backward differences (Engle, Skeel, & Drees, 2005).
     */
    int orderk;        /* related to order of approximation */
    double *sham;      /* stores shadow Hamiltonian approximations */

    MD_Dvec **dif_q;   /* backward difference table for position */
    MD_Dvec **dif_p;   /* backward difference table for momentum */
    double *dif_beta;  /* backward difference table for beta */
    double **Atab;
    double **coef;
    MD_Dvec *old_f;    /* stores previous forces */
    /* end shadow Hamiltonian */

    /*
     * Variables for weak-coupling to temperature bath.
     */
    double dt_tau;     /* timestep / relaxation_time */
    double twoke_bath;

    /*
     * Variables for Nose-Hoover.
     */
    double twice_desired_ke;  /* twice desired kinetic energy,
                                 internal units of AMU A^2 / fs^2 */
    double nh_qmass;          /* thermostat "mass" for extended Lagrangian */
    double dt_over_qmass;     /* timestep / Qmass */
    double xi;                /* thermodynamic friction coefficient */
    double eta;               /* log s */
    double twoke;             /* twice kinetic energy */
    
    /*
     * Variables for thermalized Drude oscillators.
     *
     * Atoms are transformed into
     * com (centers of mass) and dbond (Drude bonds)
     * which are propagated using a dual thermostat,
     * warm for com, cold for dbond.
     *
     * natoms = 2*num_drudeBond + num_nonPolarized
     *
     * num_dbond = num_drudeBond
     * num_com = num_drudeBond + num_nonPolarized
     */

    /* this is transformed system state */
    int32 num_com, num_dbond;
    MD_Dvec *pos_com, *pos_dbond;
    MD_Dvec *vel_com, *vel_dbond;
    MD_Dvec *f_com, *f_dbond;
    double *mass_com, *mass_dbond;
    double *scal_inv_mass_com, *scal_inv_mass_dbond;

    /* this is for dual thermostat */
    double eta_com, eta_dbond;
    double xi_com, xi_dbond;
    double qmass_com, qmass_dbond;
    double twice_desired_ke_com;
    double twice_desired_ke_dbond;
    double dt_over_qmass_com;
    double dt_over_qmass_dbond;
    double twoke_com, twoke_dbond;
    double tempkonst_com, tempkonst_dbond;

    /* this is the map that transforms between atoms and (com, dbond) */
    DrudeBond *drudeBond;    /* the Drude bond topology */
    int32 num_drudeBond;
    int32 *nonPolarized;     /* the remaining non-polarized atoms */
    int32 num_nonPolarized;
    int32 is_drude;

    /*
     * Variables for CG minimization.
     */
    double u;
    double fdf;
    MD_Dvec *starting_position;
    MD_Dvec *search_direction;
    /* end CG minimization */

    /*
     * Variables for SETTLE
     */
    MD_Dvec *refpos;  /* have to retain reference positions */
    double o_mass;  /* mass of oxygen */
    double h_mass;  /* mass of hydrogen */
    double ra, rb, rc;  /* canonical positions of water atoms */
  };


  /* for error handling */
  int step_error(Step *, const char *fmt, ...);

  /* for setting up velocities */
  int step_set_random_velocities(Step *);
  void step_remove_com_motion(Step *);

  /* callback to provide simulation results */
  int step_results(Step *, int32 delta_stepnum);
  void step_results_twoke(double *e_sum,
      const MD_Dvec *v, const MD_Atom *atom, int32 natoms);
  void step_results_linmo(MD_Dvec *linmo,
      const MD_Dvec *v, const MD_Atom *atom, int32 natoms);
  void step_results_angmo(MD_Dvec *angmo,
      const MD_Dvec *r, const MD_Dvec *v, const MD_Atom *atom, int32 natoms);
  void step_results_kvirial(double u_kv[VIRIAL_UPPER_LEN],
      const MD_Dvec *v, const MD_Atom *atom, int32 natoms);

  /* callback to compute force */
  int step_force(Step *);

  /* callback to send text output */
  int step_output(Step *, const char *fmt, ...);

  /* integration methods */
  int step_compute_verlet(Step *, int32 numsteps);
  int step_init_verlet(Step *);
  void step_done_verlet(Step *);

  int step_compute_shadow(Step *, int32 numsteps);
  int step_init_shadow(Step *);
  void step_done_shadow(Step *);

  int step_compute_tempbath(Step *, int32 numsteps);
  int step_init_tempbath(Step *);
  void step_done_tempbath(Step *);

  int step_compute_nosehoover_explicit(Step *, int32 numsteps);
  int step_init_nosehoover_explicit(Step *);
  void step_done_nosehoover_explicit(Step *);
  
  int step_compute_drude_thermal(Step *, int32 numsteps);
  int step_init_drude_thermal(Step *);
  void step_done_drude_thermal(Step *);
  
  int step_compute_drude_roux(Step *, int32 numsteps);
  int step_init_drude_roux(Step *);
  void step_done_drude_roux(Step *);
  
  int step_compute_drude_chen(Step *, int32 numsteps);
  int step_init_drude_chen(Step *);
  void step_done_drude_chen(Step *);
  
  int step_compute_cgmin(Step *, int32 numsteps);
  int step_init_cgmin(Step *);
  void step_done_cgmin(Step *);


  /* for SETTLE */
  int step_setup_settle(Step *s);
  void step_cleanup_settle(Step *s);
  int step_settle_startup(Step *s);
  int step_settle_prep(Step *s);
  int step_settle1(Step *s, const double dt);
  int step_settle2(Step *s, const double dt);


#ifdef __cplusplus
}
#endif

#endif /* STEPDEFN_H */

---