ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/SC.cpp
(Generate patch)

Comparing branches/development/src/nonbonded/SC.cpp (file contents):
Revision 1489 by gezelter, Tue Aug 10 18:34:59 2010 UTC vs.
Revision 1549 by gezelter, Wed Apr 27 18:38:15 2011 UTC

# Line 49 | Line 49 | namespace OpenMD {
49  
50   namespace OpenMD {
51  
52 –  bool SC::initialized_ = false;
53 –  RealType SC::scRcut_ = 0.0;
54 –  int SC::np_ = 3000;
55 –  ForceField* SC::forceField_ = NULL;
56 –  map<int, AtomType*> SC::SClist;
57 –  map<AtomType*, SCAtomData> SC::SCMap;
58 –  map<pair<AtomType*, AtomType*>, SCInteractionData> SC::MixingMap;
59 –  
60 –  SC* SC::_instance = NULL;
52  
53 <  SC* SC::Instance() {
54 <    if (!_instance) {
64 <      _instance = new SC();
65 <    }
66 <    return _instance;
67 <  }
53 >  SC::SC() : name_("SC"), initialized_(false), forceField_(NULL),
54 >             scRcut_(0.0), np_(3000) {}
55    
56    SCParam SC::getSCParam(AtomType* atomType) {
57      
# Line 315 | Line 302 | namespace OpenMD {
302      return;
303    }
304  
305 <  void SC::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
319 <                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
305 >  void SC::calcDensity(InteractionData &idat) {
306      
307      if (!initialized_) initialize();
308      
309 <    SCInteractionData mixer = MixingMap[make_pair(at1, at2)];
309 >    SCInteractionData mixer = MixingMap[idat.atypes];
310  
311 <    rho_i_at_j = mixer.phi->getValueAt(rij);
326 <    rho_j_at_i = rho_i_at_j;
311 >    RealType rcij = mixer.rCut;
312  
313 +    if (idat.rij < rcij) {
314 +      idat.rho_i_at_j = mixer.phi->getValueAt(idat.rij);
315 +      idat.rho_j_at_i = idat.rho_i_at_j;
316 +    } else {
317 +      idat.rho_i_at_j = 0.0;
318 +      idat.rho_j_at_i = 0.0;
319 +    }
320 +
321      return;
322    }
323  
324 <  void SC::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
332 <                          RealType &dfrhodrho) {
324 >  void SC::calcFunctional(SelfData &sdat) {
325  
326      if (!initialized_) initialize();
327  
328 <    SCAtomData data1 = SCMap[at1];
328 >    SCAtomData data1 = SCMap[sdat.atype];
329      
330 <    frho = - data1.c * data1.epsilon * sqrt(rho);
331 <    dfrhodrho = 0.5 * frho / rho;
330 >    sdat.frho = - data1.c * data1.epsilon * sqrt(sdat.rho);
331 >    sdat.dfrhodrho = 0.5 * sdat.frho / sdat.rho;
332      
333      return;
334    }
335 <
335 >  
336  
337 <  void SC::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
346 <                     RealType rij, RealType r2, RealType sw,
347 <                     RealType &vpair, RealType &pot, Vector3d &f1,
348 <                     RealType rho_i, RealType rho_j,
349 <                     RealType dfrhodrho_i, RealType dfrhodrho_j,
350 <                     RealType &fshift_i, RealType &fshift_j) {
337 >  void SC::calcForce(InteractionData &idat) {
338      
339      if (!initialized_) initialize();
340      
341 <    SCAtomData data1 = SCMap[at1];
342 <    SCAtomData data2 = SCMap[at1];
341 >    SCAtomData data1 = SCMap[idat.atypes.first];
342 >    SCAtomData data2 = SCMap[idat.atypes.second];
343  
344 <    SCInteractionData mixer = MixingMap[make_pair(at1, at2)];
344 >    SCInteractionData mixer = MixingMap[idat.atypes];
345  
346      RealType rcij = mixer.rCut;
360 –    RealType vcij = mixer.vCut;
347  
348 <    pair<RealType, RealType> res;
349 <
350 <    res = mixer.phi->getValueAndDerivativeAt(rij);
351 <    RealType rhtmp = res.first;
352 <    RealType drhodr = res.second;
353 <
354 <    res = mixer.V->getValueAndDerivativeAt(rij);
355 <    RealType vptmp = res.first;
356 <    RealType dvpdr = res.second;
357 <
358 <    RealType pot_temp = vptmp - vcij;
359 <    vpair += pot_temp;
360 <
361 <    RealType dudr = drhodr * (dfrhodrho_i + dfrhodrho_j) + dvpdr;
362 <
363 <    f1 += d * dudr / rij;
348 >    if (idat.rij < rcij) {
349 >      RealType vcij = mixer.vCut;
350 >      
351 >      pair<RealType, RealType> res;
352 >      
353 >      res = mixer.phi->getValueAndDerivativeAt(idat.rij);
354 >      RealType rhtmp = res.first;
355 >      RealType drhodr = res.second;
356 >      
357 >      res = mixer.V->getValueAndDerivativeAt(idat.rij);
358 >      RealType vptmp = res.first;
359 >      RealType dvpdr = res.second;
360 >      
361 >      RealType pot_temp = vptmp - vcij;
362 >      idat.vpair += pot_temp;
363 >      
364 >      RealType dudr = drhodr * (idat.dfrho1 + idat.dfrho2) + dvpdr;
365 >      
366 >      idat.f1 += idat.d * dudr / idat.rij;
367          
368 <    // particle_pot is the difference between the full potential
369 <    // and the full potential without the presence of a particular
370 <    // particle (atom1).
371 <    //
372 <    // This reduces the density at other particle locations, so
373 <    // we need to recompute the density at atom2 assuming atom1
374 <    // didn't contribute.  This then requires recomputing the
375 <    // density functional for atom2 as well.
376 <    //
377 <    // Most of the particle_pot heavy lifting comes from the
378 <    // pair interaction, and will be handled by vpair.
379 <    
380 <    fshift_i = - data1.c * data1.epsilon * sqrt(rho_i - rhtmp);
381 <    fshift_j = - data2.c * data2.epsilon * sqrt(rho_j - rhtmp);
382 <    
383 <    pot += pot_temp;
384 <
368 >      // particle_pot is the difference between the full potential
369 >      // and the full potential without the presence of a particular
370 >      // particle (atom1).
371 >      //
372 >      // This reduces the density at other particle locations, so
373 >      // we need to recompute the density at atom2 assuming atom1
374 >      // didn't contribute.  This then requires recomputing the
375 >      // density functional for atom2 as well.
376 >      //
377 >      // Most of the particle_pot heavy lifting comes from the
378 >      // pair interaction, and will be handled by vpair.
379 >      
380 >      idat.fshift1 = - data1.c * data1.epsilon * sqrt(idat.rho1 - rhtmp);
381 >      idat.fshift2 = - data2.c * data2.epsilon * sqrt(idat.rho2 - rhtmp);
382 >      
383 >      idat.pot[3] += pot_temp;
384 >    }
385 >      
386      return;    
387    }
388  
389 +  RealType SC::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
390 +    if (!initialized_) initialize();  
391  
392 <  void SC::calc_sc_prepair_rho(int *atid1, int *atid2, RealType *rij,
393 <                               RealType* rho_i_at_j, RealType* rho_j_at_i){
394 <    
395 <    if (!initialized_) initialize();
396 <    
397 <    AtomType* atype1 = SClist[*atid1];
398 <    AtomType* atype2 = SClist[*atid2];
399 <    
408 <    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
409 <    
410 <    return;    
392 >    map<pair<AtomType*, AtomType*>, SCInteractionData>::iterator it;
393 >    it = MixingMap.find(atypes);
394 >    if (it == MixingMap.end())
395 >      return 0.0;
396 >    else  {
397 >      SCInteractionData mixer = (*it).second;
398 >      return mixer.rCut;
399 >    }
400    }
412 –  
413 –  void SC::calc_sc_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
414 –                                 RealType *dfrhodrho) {
415 –    
416 –    if (!initialized_) initialize();
417 –    
418 –    AtomType* atype1 = SClist[*atid1];  
419 –
420 –    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
421 –    
422 –    return;    
423 –  }
424 –  
425 –  RealType SC::getSCcut(int *atid1) {
426 –
427 –    if (!initialized_) initialize();
428 –    
429 –    AtomType* atype1 = SClist[*atid1];  
430 –      
431 –    return 2.0 * getAlpha(atype1);
432 –  }
433 –
434 –  void SC::do_sc_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
435 –                        RealType *r2, RealType *sw, RealType *vpair,
436 –                        RealType *pot, RealType *f1, RealType *rho1,
437 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
438 –                        RealType *fshift1, RealType *fshift2) {
439 –
440 –    if (!initialized_) initialize();
441 –    
442 –    AtomType* atype1 = SClist[*atid1];
443 –    AtomType* atype2 = SClist[*atid2];
444 –    
445 –    Vector3d disp(d[0], d[1], d[2]);
446 –    Vector3d frc(f1[0], f1[1], f1[2]);
447 –    
448 –    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
449 –              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
450 –      
451 –    f1[0] = frc.x();
452 –    f1[1] = frc.y();
453 –    f1[2] = frc.z();
454 –
455 –    return;    
456 –  }
457 –  
458 –  void SC::setCutoffSC(RealType *thisRcut) {
459 –    scRcut_ = *thisRcut;
460 –  }
401   }
462 –
463 – extern "C" {
464 –  
465 – #define fortranCalcDensity FC_FUNC(calc_sc_prepair_rho, CALC_SC_PREPAIR_RHO)
466 – #define fortranCalcFunctional FC_FUNC(calc_sc_preforce_frho, CALC_SC_PREFORCE_FRHO)
467 – #define fortranCalcForce FC_FUNC(do_sc_pair, DO_SC_PAIR)
468 – #define fortranSetCutoffSC FC_FUNC(setcutoffsc, SETCUTOFFSC)
469 – #define fortranGetSCcut FC_FUNC(getsccut, GETSCCUT)
470 –  
471 –  
472 –  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij,
473 –                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
474 –    
475 –    return OpenMD::SC::Instance()->calc_sc_prepair_rho(atid1, atid2, rij,
476 –                                                       rho_i_at_j,  
477 –                                                       rho_j_at_i);
478 –  }
479 –  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
480 –                             RealType *dfrhodrho) {  
481 –    
482 –    return OpenMD::SC::Instance()->calc_sc_preforce_Frho(atid1, rho, frho,
483 –                                                         dfrhodrho);
484 –    
485 –  }
486 –  void fortranSetCutoffSC(RealType *rcut) {
487 –    return OpenMD::SC::Instance()->setCutoffSC(rcut);
488 –  }
489 –  void fortranCalcForce(int *atid1, int *atid2, RealType *d, RealType *rij,
490 –                        RealType *r2, RealType *sw, RealType *vpair,
491 –                        RealType *pot, RealType *f1, RealType *rho1,
492 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
493 –                        RealType *fshift1, RealType *fshift2){
494 –    
495 –    return OpenMD::SC::Instance()->do_sc_pair(atid1, atid2, d, rij,
496 –                                              r2, sw, vpair,
497 –                                              pot, f1, rho1,
498 –                                              rho2, dfrho1, dfrho2,
499 –                                              fshift1,  fshift2);
500 –  }
501 –  RealType fortranGetSCcut(int* atid) {
502 –    return OpenMD::SC::Instance()->getSCcut(atid);
503 –  }
504 –  
505 – }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines