1 #include "ChargeRegulator.h" 2 #include "PoissonSolver.h" 3 #include "LammpsInterface.h" 4 #include "ATC_Coupling.h" 5 #include "ATC_Error.h" 6 7 #include "Function.h" 8 #include "PrescribedDataManager.h" 9 10 #include <sstream> 11 #include <string> 12 #include <vector> 13 #include <utility> 14 #include <set> 15 16 using ATC_Utility::to_string; 17 using std::stringstream; 18 using std::map; 19 using std::vector; 20 using std::set; 21 using std::pair; 22 using std::string; 23 24 25 namespace ATC { 26 27 //======================================================== 28 // Class ChargeRegulator 29 //======================================================== 30 31 //-------------------------------------------------------- 32 // Constructor 33 //-------------------------------------------------------- ChargeRegulator(ATC_Coupling * atc)34 ChargeRegulator::ChargeRegulator(ATC_Coupling * atc) : 35 AtomicRegulator(atc) 36 { 37 // do nothing 38 } 39 //-------------------------------------------------------- 40 // Destructor 41 //-------------------------------------------------------- ~ChargeRegulator()42 ChargeRegulator::~ChargeRegulator() 43 { 44 map<string,ChargeRegulatorMethod *>::iterator it; 45 for (it = regulators_.begin(); it != regulators_.end(); it++) { 46 if (it->second) delete it->second; 47 } 48 } 49 50 51 //-------------------------------------------------------- 52 // modify: 53 // parses and adjusts charge regulator state based on 54 // user input, in the style of LAMMPS user input 55 //-------------------------------------------------------- modify(int narg,char ** arg)56 bool ChargeRegulator::modify(int narg, char **arg) 57 { 58 bool foundMatch = false; 59 return foundMatch; 60 } 61 62 //-------------------------------------------------------- 63 // construct methods 64 //-------------------------------------------------------- construct_methods()65 void ChargeRegulator::construct_methods() 66 { 67 AtomicRegulator::construct_methods(); 68 69 if (atc_->reset_methods()) { 70 // eliminate existing methods 71 delete_method(); 72 // consruct new ones 73 map<string, ChargeRegulatorParameters>::iterator itr; 74 for (itr = parameters_.begin(); 75 itr != parameters_.end(); itr++) { 76 string tag = itr->first; 77 if (regulators_.find(tag) != regulators_.end()) delete regulators_[tag]; 78 ChargeRegulatorParameters & p = itr->second; 79 LammpsInterface * lammpsInterface = LammpsInterface::instance(); 80 p.groupBit = lammpsInterface->group_bit(tag); 81 if (! p.groupBit) 82 throw ATC_Error("ChargeRegulator::initialize group not found"); 83 switch (p.method) { 84 case NONE: { 85 regulators_[tag] = new ChargeRegulatorMethod(this,p); 86 break; 87 } 88 case FEEDBACK: { 89 regulators_[tag] = new ChargeRegulatorMethodFeedback(this,p); 90 break; 91 } 92 case IMAGE_CHARGE: { 93 regulators_[tag] = new ChargeRegulatorMethodImageCharge(this,p); 94 break; 95 } 96 case EFFECTIVE_CHARGE: { 97 regulators_[tag] = new ChargeRegulatorMethodEffectiveCharge(this,p); 98 break; 99 } 100 default: 101 throw ATC_Error("ChargeRegulator::construct_method unknown charge regulator type"); 102 } 103 } 104 } 105 } 106 107 //-------------------------------------------------------- 108 // initialize: 109 //-------------------------------------------------------- initialize()110 void ChargeRegulator::initialize() 111 { 112 113 114 map<string, ChargeRegulatorMethod *>::iterator itr; 115 for (itr = regulators_.begin(); 116 itr != regulators_.end(); itr++) { itr->second->initialize(); } 117 118 atc_->set_boundary_integration_type(boundaryIntegrationType_); 119 AtomicRegulator::reset_nlocal(); 120 AtomicRegulator::delete_unused_data(); 121 needReset_ = false; 122 } 123 124 //-------------------------------------------------------- 125 // apply pre force 126 //-------------------------------------------------------- apply_pre_force(double dt)127 void ChargeRegulator::apply_pre_force(double dt) 128 { 129 map<string, ChargeRegulatorMethod *>::iterator itr; 130 for (itr = regulators_.begin(); 131 itr != regulators_.end(); itr++) { itr->second->apply_pre_force(dt);} 132 } 133 //-------------------------------------------------------- 134 // apply post force 135 //-------------------------------------------------------- apply_post_force(double dt)136 void ChargeRegulator::apply_post_force(double dt) 137 { 138 map<string, ChargeRegulatorMethod *>::iterator itr; 139 for (itr = regulators_.begin(); 140 itr != regulators_.end(); itr++) { itr->second->apply_post_force(dt);} 141 } 142 //-------------------------------------------------------- 143 // output 144 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)145 void ChargeRegulator::output(OUTPUT_LIST & outputData) 146 { 147 map<string, ChargeRegulatorMethod *>::iterator itr; 148 for (itr = regulators_.begin(); 149 itr != regulators_.end(); itr++) { itr->second->output(outputData);} 150 } 151 152 //======================================================== 153 // Class ChargeRegulatorMethod 154 //======================================================== 155 156 //-------------------------------------------------------- 157 // Constructor 158 // Grab references to ATC and ChargeRegulator 159 //-------------------------------------------------------- ChargeRegulatorMethod(ChargeRegulator * chargeRegulator,ChargeRegulator::ChargeRegulatorParameters & p)160 ChargeRegulatorMethod::ChargeRegulatorMethod 161 (ChargeRegulator *chargeRegulator, 162 ChargeRegulator::ChargeRegulatorParameters & p) 163 : RegulatorShapeFunction(chargeRegulator), 164 chargeRegulator_(chargeRegulator), 165 lammpsInterface_(LammpsInterface::instance()), 166 rC_(0), rCsq_(0), 167 targetValue_(NULL), 168 targetPhi_(p.value), 169 surface_(p.faceset), 170 atomGroupBit_(p.groupBit), 171 boundary_(false), 172 depth_(p.depth), 173 surfaceType_(p.surfaceType), 174 permittivity_(p.permittivity), 175 initialized_(false) 176 { 177 const FE_Mesh * feMesh = atc_->fe_engine()->fe_mesh(); 178 feMesh->faceset_to_nodeset(surface_,nodes_); 179 // assume flat get normal and primary coord 180 PAIR face = *(surface_.begin()); 181 normal_.reset(nsd_); 182 feMesh->face_normal(face,0,normal_); 183 DENS_MAT faceCoords; 184 feMesh->face_coordinates(face,faceCoords); 185 point_.reset(nsd_); 186 for (int i=0; i < nsd_; i++) { point_(i) = faceCoords(i,0); } 187 #ifdef ATC_VERBOSE 188 stringstream ss; ss << "point: (" << point_(0) << "," << point_(1) << "," << point_(2) << ") normal: (" << normal_(0) << "," << normal_(1) << "," << normal_(2) << ") depth: " << depth_; 189 lammpsInterface_->print_msg_once(ss.str()); 190 #endif 191 sum_.reset(nsd_); 192 } 193 194 //-------------------------------------------------------- 195 // Initialize 196 197 //-------------------------------------------------------- 198 199 200 // nomenclature might be a bit backwark: control --> nodes that exert the control, & influence --> atoms that feel the influence initialize(void)201 void ChargeRegulatorMethod::initialize(void) 202 { 203 interscaleManager_ = &(atc_->interscale_manager()); 204 205 poissonSolver_ =chargeRegulator_->poisson_solver(); 206 if (! poissonSolver_) throw ATC_Error("need a poisson solver to initialize charge regulator"); 207 208 // atomic vectors 209 // nodal information 210 nNodes_ = atc_->num_nodes(); 211 // constants 212 rC_ = lammpsInterface_->pair_cutoff(); 213 rCsq_ = rC_*rC_; 214 qV2e_ = lammpsInterface_->qv2e(); 215 qqrd2e_ = lammpsInterface_->qqrd2e(); 216 217 // note derived method set intialized to true 218 } 219 220 nlocal()221 int ChargeRegulatorMethod::nlocal() { return atc_->nlocal(); } 222 set_greens_functions(void)223 void ChargeRegulatorMethod::set_greens_functions(void) 224 { 225 // set up Green's function per node 226 for (int i = 0; i < nNodes_; i++) { 227 set<int> localNodes; 228 for (int j = 0; j < nNodes_; j++) 229 localNodes.insert(j); 230 // call Poisson solver to get Green's function for node i 231 DENS_VEC globalGreensFunction; 232 poissonSolver_->greens_function(i,globalGreensFunction); 233 // store green's functions as sparse vectors only on local nodes 234 set<int>::const_iterator thisNode; 235 SparseVector<double> sparseGreensFunction(nNodes_); 236 for (thisNode = localNodes.begin(); thisNode != localNodes.end(); thisNode++) 237 sparseGreensFunction(*thisNode) = globalGreensFunction(*thisNode); 238 greensFunctions_.push_back(sparseGreensFunction); 239 } 240 } 241 //-------------------------------------------------------- 242 // output 243 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)244 void ChargeRegulatorMethod::output(OUTPUT_LIST & outputData) 245 { 246 //vector<double> localSum(sum_.size()); 247 //lammpsInteface_->allsum(localSum.pointer,sum_.pointer,sum_.size()); 248 DENS_VEC localSum(sum_.size()); 249 lammpsInterface_->allsum(localSum.ptr(),sum_.ptr(),sum_.size()); 250 for (int i = 0; i < sum_.size(); i++) { 251 string name = "charge_regulator_influence_"+to_string(i); 252 // atc_->fe_engine()->add_global(name,sum_[i]); 253 } 254 } 255 256 //======================================================== 257 // Class ChargeRegulatorMethodFeedback 258 //======================================================== 259 //-------------------------------------------------------- 260 // Constructor 261 //-------------------------------------------------------- ChargeRegulatorMethodFeedback(ChargeRegulator * chargeRegulator,ChargeRegulator::ChargeRegulatorParameters & p)262 ChargeRegulatorMethodFeedback::ChargeRegulatorMethodFeedback 263 (ChargeRegulator *chargeRegulator, 264 ChargeRegulator::ChargeRegulatorParameters & p) 265 : ChargeRegulatorMethod (chargeRegulator, p), 266 controlNodes_(nodes_), 267 influenceGroupBit_(p.groupBit) 268 { 269 nControlNodes_ = controlNodes_.size(); 270 sum_.resize(1); 271 } 272 //-------------------------------------------------------- 273 // Initialize 274 //-------------------------------------------------------- initialize(void)275 void ChargeRegulatorMethodFeedback::initialize(void) 276 { 277 ChargeRegulatorMethod::initialize(); 278 if (surfaceType_ != ChargeRegulator::CONDUCTOR) 279 throw ATC_Error("currently charge feedback can only mimic a conductor"); 280 set_influence(); 281 set_influence_matrix(); 282 initialized_ = true; 283 } 284 //-------------------------------------------------------- 285 // Initialize 286 //-------------------------------------------------------- construct_transfers(void)287 void ChargeRegulatorMethodFeedback::construct_transfers(void) 288 { 289 ChargeRegulatorMethod::construct_transfers(); 290 291 InterscaleManager & interscaleManager((atomicRegulator_->atc_transfer())->interscale_manager()); 292 PerAtomSparseMatrix<double> * atomShapeFunctions = interscaleManager.per_atom_sparse_matrix("InterpolantGhost"); 293 if (!atomShapeFunctions) { 294 atomShapeFunctions = new PerAtomShapeFunction(atomicRegulator_->atc_transfer(), 295 interscaleManager.per_atom_quantity("AtomicGhostCoarseGrainingPositions"), 296 interscaleManager.per_atom_int_quantity("AtomGhostElement"), 297 GHOST); 298 interscaleManager.add_per_atom_sparse_matrix(atomShapeFunctions,"InterpolantGhost"); 299 } 300 } 301 //-------------------------------------------------------- 302 // find measurement atoms and nodes 303 //-------------------------------------------------------- set_influence(void)304 void ChargeRegulatorMethodFeedback::set_influence(void) 305 { 306 307 // get nodes that overlap influence atoms & compact list of influence atoms 308 boundary_ = 309 atc_->nodal_influence(influenceGroupBit_,influenceNodes_,influenceAtoms_); 310 nInfluenceAtoms_ = influenceAtoms_.size(); // local 311 nInfluenceNodes_ = influenceNodes_.size(); // global 312 stringstream ss; ss << "control nodes: " << nControlNodes_ << " influence nodes: " << nInfluenceNodes_ << " local influence atoms: " << nInfluenceAtoms_ ; 313 lammpsInterface_->print_msg(ss.str()); 314 if (nInfluenceNodes_ == 0) throw ATC_Error("no influence nodes"); 315 316 const Array<int> & map = (boundary_) ? atc_->ghost_to_atom_map() : atc_->internal_to_atom_map(); 317 for (set<int>::const_iterator itr = influenceAtoms_.begin(); itr != influenceAtoms_.end(); itr++) { 318 influenceAtomsIds_.insert(map(*itr)); 319 } 320 } 321 //-------------------------------------------------------- 322 // constuct a Green's submatrix 323 //-------------------------------------------------------- set_influence_matrix(void)324 void ChargeRegulatorMethodFeedback::set_influence_matrix(void) 325 { 326 // construct control-influence matrix bar{G}^-1: ds{p} = G{p,m}^-1 dphi{m} 327 328 329 // 330 if (nInfluenceNodes_ < nControlNodes_) throw ATC_Error(" least square not implmented "); 331 if (nInfluenceNodes_ > nControlNodes_) throw ATC_Error(" solve not possible "); 332 DENS_MAT G(nInfluenceNodes_,nControlNodes_); 333 DENS_VEC G_I; 334 set<int>::const_iterator itr,itr2,itr3; 335 const Array<int> & nmap = atc_->fe_engine()->fe_mesh()->global_to_unique_map(); 336 int i = 0; 337 for (itr = influenceNodes_.begin(); itr != influenceNodes_.end(); itr++) { 338 poissonSolver_->greens_function(*itr, G_I); 339 int j = 0; 340 for (itr2 = controlNodes_.begin(); itr2 != controlNodes_.end(); itr2++) { 341 int jnode = nmap(*itr2); 342 G(i,j++) = G_I(jnode); 343 } 344 i++; 345 } 346 invG_ = inv(G); 347 348 // construct the prolong-restrict projector N N^T for influence nodes only 349 350 InterscaleManager & interscaleManager(atc_->interscale_manager()); 351 const SPAR_MAT & N_Ia = (boundary_) ? 352 (interscaleManager.per_atom_sparse_matrix("InterpolantGhost"))->quantity(): 353 (interscaleManager.per_atom_sparse_matrix("Interpolant"))->quantity(); 354 NT_.reset(nInfluenceAtoms_,nInfluenceNodes_); 355 DENS_MAT NNT(nInfluenceNodes_,nInfluenceNodes_); 356 int k = 0; 357 for (itr3 = influenceAtoms_.begin(); itr3 != influenceAtoms_.end(); itr3++) { 358 int katom = *itr3; 359 int i = 0; 360 for (itr = influenceNodes_.begin(); itr != influenceNodes_.end(); itr++) { 361 int Inode = *itr; 362 int j = 0; 363 NT_(k,i) = N_Ia(katom,Inode); 364 for (itr2 = influenceNodes_.begin(); itr2 != influenceNodes_.end(); itr2++) { 365 int Jnode = *itr2; 366 NNT(i,j++) += N_Ia(katom,Inode)*N_Ia(katom,Jnode); 367 } 368 i++; 369 } 370 k++; 371 } 372 // swap contributions across processors 373 DENS_MAT localNNT = NNT; 374 int count = NNT.nRows()*NNT.nCols(); 375 lammpsInterface_->allsum(localNNT.ptr(),NNT.ptr(),count); 376 invNNT_ = inv(NNT); 377 378 // total influence matrix 379 if (nInfluenceAtoms_ > 0) { NTinvNNTinvG_ = NT_*invNNT_*invG_; } 380 381 } 382 383 //-------------------------------------------------------- 384 // change potential/charge pre-force calculation 385 //-------------------------------------------------------- apply_pre_force(double dt)386 void ChargeRegulatorMethodFeedback::apply_pre_force(double dt) 387 { 388 389 sum_ = 0; 390 if (nInfluenceAtoms_ == 0) return; // nothing to do 391 apply_feedback_charges(); 392 } 393 //-------------------------------------------------------- 394 // apply feedback charges to atoms 395 //-------------------------------------------------------- apply_feedback_charges()396 void ChargeRegulatorMethodFeedback::apply_feedback_charges() 397 { 398 double * q = lammpsInterface_->atom_charge(); 399 // calculate error in potential on the control nodes 400 401 const DENS_MAT & phiField = (atc_->field(ELECTRIC_POTENTIAL)).quantity(); 402 DENS_MAT dphi(nControlNodes_,1); 403 int i = 0; 404 set<int>::const_iterator itr; 405 for (itr = controlNodes_.begin(); itr != controlNodes_.end(); itr++) { 406 dphi(i++,0) = targetPhi_ - phiField(*itr,0); 407 } 408 409 // construct the atomic charges consistent with the correction 410 DENS_MAT dq = NTinvNNTinvG_*dphi; 411 i = 0; 412 for (itr = influenceAtomsIds_.begin(); itr != influenceAtomsIds_.end(); itr++) { 413 sum_(0) += dq(i,0); 414 q[*itr] += dq(i++,0); 415 } 416 417 (interscaleManager_->fundamental_atom_quantity(LammpsInterface::ATOM_CHARGE))->force_reset(); 418 (interscaleManager_->fundamental_atom_quantity(LammpsInterface::ATOM_CHARGE, GHOST))->force_reset(); 419 } 420 421 //======================================================== 422 // Class ChargeRegulatorMethodImageCharge 423 //======================================================== 424 //-------------------------------------------------------- 425 // Constructor 426 //-------------------------------------------------------- ChargeRegulatorMethodImageCharge(ChargeRegulator * chargeRegulator,ChargeRegulator::ChargeRegulatorParameters & p)427 ChargeRegulatorMethodImageCharge::ChargeRegulatorMethodImageCharge 428 (ChargeRegulator *chargeRegulator, 429 ChargeRegulator::ChargeRegulatorParameters & p) 430 : ChargeRegulatorMethod (chargeRegulator, p), 431 imageNodes_(nodes_) 432 { 433 } 434 //-------------------------------------------------------- 435 // Initialize 436 //-------------------------------------------------------- initialize(void)437 void ChargeRegulatorMethodImageCharge::initialize(void) 438 { 439 ChargeRegulatorMethod::initialize(); 440 if (surfaceType_ != ChargeRegulator::DIELECTRIC) throw ATC_Error("currently image charge can only mimic a dielectric"); 441 double eps1 = permittivity_;// dielectric 442 double eps2 = lammpsInterface_->dielectric();// ambient 443 permittivityRatio_ = (eps2-eps1)/(eps2+eps1); 444 #ifdef ATC_VERBOSE 445 stringstream ss; ss << "permittivity ratio: " << permittivityRatio_; 446 lammpsInterface_->print_msg_once(ss.str()); 447 #endif 448 set_greens_functions(); 449 450 /////////////////////////////////////////////////////////////// 451 /////////////////////////////////////////////////////////////// 452 initialized_ = true; 453 } 454 455 //-------------------------------------------------------- 456 // change potential/charge post-force calculation 457 //-------------------------------------------------------- apply_post_force(double dt)458 void ChargeRegulatorMethodImageCharge::apply_post_force(double dt) 459 { 460 sum_ = 0; 461 apply_local_forces(); 462 463 //correct_forces(); 464 } 465 466 //-------------------------------------------------------- 467 // apply local coulomb forces 468 // -- due to image charges 469 //-------------------------------------------------------- apply_local_forces()470 void ChargeRegulatorMethodImageCharge::apply_local_forces() 471 { 472 473 int inum = lammpsInterface_->neighbor_list_inum(); 474 int * ilist = lammpsInterface_->neighbor_list_ilist(); 475 int * numneigh = lammpsInterface_->neighbor_list_numneigh(); 476 int ** firstneigh = lammpsInterface_->neighbor_list_firstneigh(); 477 478 const int *mask = lammpsInterface_->atom_mask(); 479 ///.............................................. 480 double ** x = lammpsInterface_->xatom(); 481 double ** f = lammpsInterface_->fatom(); 482 double * q = lammpsInterface_->atom_charge(); 483 484 // loop over neighbor list 485 for (int ii = 0; ii < inum; ii++) { 486 int i = ilist[ii]; 487 double qi = q[i]; 488 if ((mask[i] & atomGroupBit_) && qi != 0.) { 489 double* fi = f[i]; 490 DENS_VEC xi(x[i],nsd_); 491 // distance to surface 492 double dn = reflect(xi); 493 // all ions near the interface/wall 494 // (a) self image 495 if (dn < rC_) { // close enough to wall to have explicit image charges 496 double factor_coul = 1; 497 double dx = 2.*dn; // distance to image charge 498 double fn = factor_coul*qi*qi*permittivityRatio_/dx; 499 fi[0] += fn*normal_[0]; 500 fi[1] += fn*normal_[1]; 501 fi[2] += fn*normal_[2]; 502 sum_ += fn*normal_; 503 // (b) neighbor images 504 int * jlist = firstneigh[i]; 505 int jnum = numneigh[i]; 506 for (int jj = 0; jj < jnum; jj++) { 507 int j = jlist[jj]; 508 // this changes j 509 double factor_coul = lammpsInterface_->coulomb_factor(j); 510 double qj = q[j]; 511 if (qj != 0.) { // all charged neighbors 512 DENS_VEC xj(x[j],nsd_); 513 dn = reflect(xj); 514 DENS_VEC dx = xi-xj; 515 double r2 = dx.norm_sq(); 516 // neighbor image j' inside cutoff from i 517 if (r2 < rCsq_) { 518 double fm = factor_coul*qi*qj*permittivityRatio_/r2; 519 fi[0] += fm*dx(0); 520 fi[1] += fm*dx(1); 521 fi[2] += fm*dx(2); 522 sum_ += fm*dx; 523 } 524 } 525 } 526 } // end i < rC if 527 } 528 } 529 // update managed data 530 (interscaleManager_->fundamental_atom_quantity(LammpsInterface::ATOM_FORCE))->force_reset(); 531 } 532 533 //-------------------------------------------------------- 534 // correct charge densities 535 // - to reflect image charges 536 //-------------------------------------------------------- correct_charge_densities()537 void ChargeRegulatorMethodImageCharge::correct_charge_densities() 538 { 539 } 540 541 542 //-------------------------------------------------------- 543 // correct_forces 544 // - due to image charge density used in short-range solution 545 //-------------------------------------------------------- correct_forces()546 void ChargeRegulatorMethodImageCharge::correct_forces() 547 { 548 } 549 550 551 //======================================================== 552 // Class ChargeRegulatorMethodEffectiveCharge 553 //======================================================== 554 //-------------------------------------------------------- 555 // Constructor 556 //-------------------------------------------------------- 557 ChargeRegulatorMethodEffectiveCharge(ChargeRegulator * chargeRegulator,ChargeRegulator::ChargeRegulatorParameters & p)558 ChargeRegulatorMethodEffectiveCharge::ChargeRegulatorMethodEffectiveCharge( 559 ChargeRegulator *chargeRegulator, 560 ChargeRegulator::ChargeRegulatorParameters & p) 561 : ChargeRegulatorMethod (chargeRegulator, p), 562 chargeDensity_(p.value), 563 useSlab_(false) 564 { 565 } 566 //-------------------------------------------------------- 567 // add_charged_surface 568 //-------------------------------------------------------- initialize()569 void ChargeRegulatorMethodEffectiveCharge::initialize( ) 570 { 571 ChargeRegulatorMethod::initialize(); 572 boundary_ = atc_->is_ghost_group(atomGroupBit_); 573 // set face sources to all point at unit function for use in integration 574 SURFACE_SOURCE faceSources; 575 map<PAIR, Array<XT_Function*> > & fs(faceSources[ELECTRIC_POTENTIAL]); 576 XT_Function * f = XT_Function_Mgr::instance()->constant_function(1.); 577 set< PAIR >::const_iterator fsItr; 578 for (fsItr = surface_.begin(); fsItr != surface_.end(); fsItr++) { 579 Array < XT_Function * > & dof = fs[*fsItr]; 580 dof.reset(1); 581 dof(0) = f; 582 } 583 584 // computed integrals of nodal shape functions on face 585 FIELDS nodalFaceWeights; 586 Array<bool> fieldMask(NUM_FIELDS); fieldMask(ELECTRIC_POTENTIAL) = true; 587 (atc_->fe_engine())->compute_fluxes(fieldMask,0.,faceSources,nodalFaceWeights); 588 const DENS_MAT & w = (nodalFaceWeights[ELECTRIC_POTENTIAL].quantity()); 589 590 // Get coordinates of each node in face set 591 for (set<int>::const_iterator n =nodes_.begin(); n != nodes_.end(); n++) { 592 DENS_VEC x = atc_->fe_engine()->fe_mesh()->nodal_coordinates(*n); 593 // compute effective charge at each node I 594 // multiply charge density by integral of N_I over face 595 double v = w(*n,0)*chargeDensity_; 596 pair<DENS_VEC,double> p(x,v); 597 nodeXFMap_[*n] = p; 598 } 599 600 // set up data structure holding charged faceset information 601 FIELDS sources; 602 double k = lammpsInterface_->coulomb_constant(); 603 string fname = "radial_power"; 604 double xtArgs[8]; 605 xtArgs[0] = 0; xtArgs[1] = 0; xtArgs[2] = 0; 606 xtArgs[3] = 1; xtArgs[4] = 1; xtArgs[5] = 1; 607 xtArgs[6] = k*chargeDensity_; 608 xtArgs[7] = -1.; 609 const DENS_MAT & s(sources[ELECTRIC_POTENTIAL].quantity()); 610 NODE_TO_XF_MAP::iterator XFitr; 611 for (XFitr = nodeXFMap_.begin(); XFitr != nodeXFMap_.end(); XFitr++) { 612 // evaluate voltage at each node I 613 // set up X_T function for integration: k*chargeDensity_/||x_I - x_s|| 614 // integral is approximated in two parts: 615 // 1) near part with all faces within r < rcrit evaluated as 2 * pi * rcrit * k sigma A/A0, A is area of this region and A0 = pi * rcrit^2, so 2 k sigma A / rcrit 616 // 2) far part evaluated using Gaussian quadrature on faceset 617 DENS_VEC x((XFitr->second).first); 618 xtArgs[0] = x(0); xtArgs[1] = x(1); xtArgs[2] = x(2); 619 f = XT_Function_Mgr::instance()->function(fname,8,xtArgs); 620 for (fsItr = surface_.begin(); fsItr != surface_.end(); fsItr++) { 621 fs[*fsItr] = f; 622 } 623 624 // perform integration to get quantities at nodes on facesets 625 // V_J' = int_S N_J k*sigma/|x_I - x_s| dS 626 (atc_->fe_engine())->compute_fluxes(fieldMask,0.,faceSources,sources); 627 628 // sum over all nodes in faceset to get total potential: 629 // V_I = sum_J VJ' 630 int node = XFitr->first; 631 nodalChargePotential_[node] = s(node,0); 632 double totalPotential = 0.; 633 for (set<int>::const_iterator n =nodes_.begin(); n != nodes_.end(); n++) { 634 totalPotential += s(*n,0); } 635 636 // assign an XT function per each node and 637 // then call the prescribed data manager and fix each node individually. 638 f = XT_Function_Mgr::instance()->constant_function(totalPotential); 639 (atc_->prescribed_data_manager())->fix_field(node,ELECTRIC_POTENTIAL,0,f); 640 } 641 initialized_ = true; 642 } 643 644 //-------------------------------------------------------- 645 // add effective forces post LAMMPS force call 646 //-------------------------------------------------------- apply_post_force(double dt)647 void ChargeRegulatorMethodEffectiveCharge::apply_post_force(double dt) 648 { 649 apply_local_forces(); 650 } 651 652 //-------------------------------------------------------- 653 // apply_charged_surfaces 654 //-------------------------------------------------------- apply_local_forces()655 void ChargeRegulatorMethodEffectiveCharge::apply_local_forces() 656 { 657 double * q = lammpsInterface_->atom_charge(); 658 _atomElectricalForce_.resize(nlocal(),nsd_); 659 660 double penalty = poissonSolver_->penalty_coefficient(); 661 if (penalty <= 0.0) throw ATC_Error("ExtrinsicModelElectrostatic::apply_charged_surfaces expecting non zero penalty"); 662 663 double dx[3]; 664 const DENS_MAT & xa((interscaleManager_->per_atom_quantity("AtomicCoarseGrainingPositions"))->quantity()); 665 666 // WORKSPACE - most are static 667 SparseVector<double> dv(nNodes_); 668 vector<SparseVector<double> > derivativeVectors; 669 derivativeVectors.reserve(nsd_); 670 const SPAR_MAT_VEC & shapeFunctionDerivatives((interscaleManager_->vector_sparse_matrix("InterpolateGradient"))->quantity()); 671 672 DenseVector<INDEX> nodeIndices; 673 DENS_VEC nodeValues; 674 675 NODE_TO_XF_MAP::const_iterator inode; 676 for (inode = nodeXFMap_.begin(); inode != nodeXFMap_.end(); inode++) { 677 678 int node = inode->first; 679 DENS_VEC xI = (inode->second).first; 680 double qI = (inode->second).second; 681 double phiI = nodalChargePotential_[node]; 682 for (int i = 0; i < nlocal(); i++) { 683 int atom = (atc_->internal_to_atom_map())(i); 684 double qa = q[atom]; 685 if (qa != 0) { 686 double dxSq = 0.; 687 for (int j = 0; j < nsd_; j++) { 688 dx[j] = xa(i,j) - xI(j); 689 dxSq += dx[j]*dx[j]; 690 } 691 if (dxSq < rCsq_) { 692 // first apply pairwise coulombic interaction 693 if (!useSlab_) { 694 double coulForce = qqrd2e_*qI*qa/(dxSq*sqrtf(dxSq)); 695 for (int j = 0; j < nsd_; j++) { 696 _atomElectricalForce_(i,j) += dx[j]*coulForce; } 697 } 698 699 // second correct for FE potential induced by BCs 700 // determine shape function derivatives at atomic location 701 // and construct sparse vectors to store derivative data 702 703 704 for (int j = 0; j < nsd_; j++) { 705 shapeFunctionDerivatives[j]->row(i,nodeValues,nodeIndices); 706 derivativeVectors.push_back(dv); 707 for (int k = 0; k < nodeIndices.size(); k++) { 708 derivativeVectors[j](nodeIndices(k)) = nodeValues(k); } 709 } 710 711 // compute greens function from charge quadrature 712 713 SparseVector<double> shortFePotential(nNodes_); 714 shortFePotential.add_scaled(greensFunctions_[node],penalty*phiI); 715 716 // compute electric field induced by charge 717 DENS_VEC efield(nsd_); 718 for (int j = 0; j < nsd_; j++) { 719 efield(j) = -.1*dot(derivativeVectors[j],shortFePotential); } 720 721 // apply correction in atomic forces 722 double c = qV2e_*qa; 723 for (int j = 0; j < nsd_; j++) { 724 if ((!useSlab_) || (j==nsd_)) { 725 _atomElectricalForce_(i,j) -= c*efield(j); 726 } 727 } 728 } 729 } 730 } 731 } 732 733 } 734 735 }; // end namespace 736