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