1 #ifndef CONCENTRATION_REGULATOR_H 2 #define CONCENTRATION_REGULATOR_H 3 4 #include <map> 5 #include <string> 6 7 #include "AtomicRegulator.h" 8 #include "LammpsInterface.h" 9 10 namespace ATC { 11 12 /** 13 * @class ConcentrationRegulator 14 * @brief Manager class for atom-continuum control of charge and potential 15 */ 16 17 class ConcentrationRegulatorMethod; 18 class ConcentrationRegulator : public AtomicRegulator { 19 20 public: 21 enum ConcentrationRegulatorType {NONE=0,TRANSITION}; 22 /** parser data */ 23 struct ConcentrationRegulatorParameters { ConcentrationRegulatorParametersConcentrationRegulatorParameters24 ConcentrationRegulatorParameters(): 25 method(NONE), 26 type(0), 27 groupbit(0), 28 transitionType(0), 29 value(0), 30 frequency(0), 31 transitionInterval(0), 32 maxEnergy(0), 33 maxExchanges(0), 34 maxAttempts(0) { }; 35 ConcentrationRegulatorType method; 36 int type; 37 int groupbit; 38 int transitionType; 39 double value; 40 int frequency; 41 int transitionInterval; 42 double maxEnergy; 43 int maxExchanges; 44 int maxAttempts; 45 std::string transitionTag; 46 ESET elemset; 47 }; 48 49 /** constructor */ 50 ConcentrationRegulator(ATC_Coupling *atc); 51 /** destructor */ 52 ~ConcentrationRegulator(); 53 /** parser/modifier */ 54 virtual bool modify(int narg, char **arg); 55 /** construct methods */ 56 virtual void construct_methods(); 57 /** pre time integration */ 58 virtual void initialize(); 59 60 //WIP_JAT need a nicer way to consistently handle sets of regulators, not sure how yet 61 // application steps 62 /** apply the regulator in the pre-predictor phase */ apply_pre_predictor(double dt,int timeStep)63 virtual void apply_pre_predictor(double dt, int timeStep){}; 64 /** apply the regulator in the mid-predictor phase */ apply_mid_predictor(double dt,int timeStep)65 virtual void apply_mid_predictor(double dt, int timeStep){}; 66 /** apply the regulator in the post-predictor phase */ apply_post_predictor(double dt,int timeStep)67 virtual void apply_post_predictor(double dt, int timeStep){}; 68 /** apply the regulator in the pre-correction phase */ apply_pre_corrector(double dt,int timeStep)69 virtual void apply_pre_corrector(double dt, int timeStep){}; 70 /** apply the regulator in the post-correction phase */ apply_post_corrector(double dt,int timeStep)71 virtual void apply_post_corrector(double dt, int timeStep){}; 72 /** compute the thermal boundary flux, must be consistent with regulator */ compute_boundary_flux(FIELDS & fields)73 virtual void compute_boundary_flux(FIELDS & fields){}; 74 /** add contributions (if any) to the finite element right-hand side */ add_to_rhs(FIELDS & rhs)75 virtual void add_to_rhs(FIELDS & rhs){}; 76 77 /** prior to exchanges */ 78 virtual void pre_force(); 79 /** prior to exchanges */ 80 virtual void pre_exchange(); 81 /** following a run */ 82 virtual void finish(); 83 /** ATC output */ 84 virtual void output(OUTPUT_LIST & outputData) const; 85 /** scalar output */ 86 virtual double compute_vector(int n) const; 87 virtual int size_vector(int s) const; 88 needs_temperature()89 bool needs_temperature() { return regulators_.size() > 0; } 90 91 protected: 92 /** registry charge regulators */ 93 std::map<std::string,ConcentrationRegulatorMethod *> regulators_; 94 std::map<std::string,ConcentrationRegulatorParameters> parameters_; 95 96 private: 97 ConcentrationRegulator(); // DO NOT define this 98 }; 99 100 /** 101 * @class ConcentrationRegulatorMethod 102 * @brief Base class for ConcentrationRegulator algorithms 103 */ 104 105 class ConcentrationRegulatorMethod : public RegulatorShapeFunction { 106 107 public: ConcentrationRegulatorMethod(ConcentrationRegulator * c)108 ConcentrationRegulatorMethod(ConcentrationRegulator *c) 109 : RegulatorShapeFunction(c) {}; ~ConcentrationRegulatorMethod()110 ~ConcentrationRegulatorMethod() {}; initialize(void)111 void initialize(void) {}; pre_force()112 virtual void pre_force() {}; pre_exchange()113 virtual void pre_exchange() {}; finish()114 virtual void finish() {}; set_weights()115 virtual void set_weights() {}; compute_vector(int n)116 virtual double compute_vector(int n) const { return 0;} output(OUTPUT_LIST & outputData)117 virtual void output(OUTPUT_LIST & outputData){}; 118 private: 119 ConcentrationRegulatorMethod(); // DO NOT define this 120 }; 121 122 /** 123 * @class ConcentrationRegulatorMethodTransition 124 * @brief GCMC + thermodynamic integration 125 */ 126 class ConcentrationRegulatorMethodTransition : public ConcentrationRegulatorMethod { 127 128 public: 129 /** constructor */ 130 ConcentrationRegulatorMethodTransition( 131 ConcentrationRegulator *concentrationRegulator, 132 ConcentrationRegulator::ConcentrationRegulatorParameters & parameters); 133 /** destructor */ ~ConcentrationRegulatorMethodTransition()134 ~ConcentrationRegulatorMethodTransition(){ 135 if (randomNumberGenerator_) delete randomNumberGenerator_; 136 } 137 /** initialize */ 138 void initialize(void); 139 /** prior to force evaluation */ 140 virtual void pre_force(); 141 /** prior to exchanges */ 142 virtual void pre_exchange(); 143 /** "thermo" output */ 144 virtual double compute_vector(int n) const; 145 protected: 146 /** set transition state: epsilon and charge */ mask(int type,int groupbit)147 int mask(int type, int groupbit) {return 0;} 148 int count(void) const; 149 int excess(void) const; 150 double energy(int id) const; 151 double energy(double * x) const; 152 double insertion_location(DENS_VEC & x) const; 153 double deletion_id(ID_PAIR & id) const ; 154 double deletion_id_consistent(ID_PAIR & id) const ; 155 double deletion_id_free(ID_PAIR & id) const ; insertion_factor(int c)156 double insertion_factor(int c) const // a ramp 157 { 158 if (c < transitionInterval_) return ((double) c)/transitionInterval_; 159 else return 1.0; 160 } 161 void transition(); 162 bool accept(double energy, double T = 0) const; 163 bool delete_atoms(int n); 164 bool insert_atoms(int n); 165 int pick_element() const; 166 void pick_coordinates(const int elem, DENS_VEC & xi, DENS_VEC& x ) const; 167 void pick_velocity(DENS_VEC & v, const double T ) const; 168 double uniform() const; 169 double normal() const; 170 171 /** pointer to conc regulator object for data */ 172 ConcentrationRegulator * concentrationRegulator_; 173 174 /** interscale manager */ 175 class InterscaleManager * interscaleManager_; 176 177 /** interscale manager */ 178 class LammpsInterface * lammpsInterface_; 179 180 /** member data */ 181 class AtomInElementSet * list_; 182 double targetConcentration_; 183 double targetCount_; 184 ESET elemset_; 185 POTENTIAL p_; 186 RNG_POINTER randomNumberGenerator_; 187 DENS_VEC volumes_; 188 double V_; 189 double q0_; 190 int controlType_; 191 int controlIndex_; 192 int transitionType_; 193 int transitionInterval_;; 194 int transitionCounter_; 195 int nInTransition_; 196 double transitionFactor_; 197 DENS_VEC epsilon0_; 198 ID_LIST deletionIds_, insertionIds_; 199 int controlMask_; 200 int frequency_; 201 double maxEnergy_; 202 int maxExchanges_; 203 int maxAttempts_; 204 int nexchanges_; 205 bool initialized_; 206 double sigma_; 207 208 void sync_random_number_generators() const; 209 mutable int _rngUniformCounter_; // for parallel consistency 210 mutable int _rngNormalCounter_; // for parallel consistency 211 212 private: 213 ConcentrationRegulatorMethodTransition(); // DO NOT define this 214 }; 215 }; 216 #endif 217