synch

jimlutsko · Sep 4, 2023 · 5635ba4 · 5635ba4
1 parent a930045
commit 5635ba4
Show file tree

Hide file tree

Showing 3 changed files with 39 additions and 0 deletions.
diff --git a/include/Species.h b/include/Species.h
@@ -80,6 +80,8 @@ class Species
   virtual void convert_to_alias_increment(DFT_Vec &dRho) const;
   virtual void convert_to_density_increment(DFT_Vec &dRho) const;
 
+  virtual void square_and_scale_with_d2rho_dx2(DFT_Vec &ff) const {throw std::runtime_error("square_and_scale_with_d2rho_dx2 not implemented");}
+
   void turn_on_display() {density_->turn_on_display();}
   void doDisplay(string &title, string &file, void *param = NULL) const { density_->doDisplay(title,file, seq_num_, param);}
 
@@ -233,6 +235,8 @@ class FMT_Species : public Species
   virtual void convert_to_alias_deriv(DFT_Vec &dF_dRho) const;
   virtual void convert_to_alias_increment(DFT_Vec &dRho) const;
   virtual void convert_to_density_increment(DFT_Vec &dRho) const;
+
+  virtual void square_and_scale_with_d2rho_dx2(DFT_Vec &ff) const;
 
   const DFT_Vec_Complex& getWEK() const { return fmt_weighted_densities[EI()].wk();}
 

diff --git a/src/DFT.cpp b/src/DFT.cpp
@@ -162,7 +162,17 @@ double DFT::calculateFreeEnergyAndDerivatives(bool onlyFex, bool H_dot_Force)
 
 	DFT_Vec &dF = s->getDF();
 	DFT_Vec ff(dF.size());
+
+	DFT_Vec ff_save(dF);
+
+	s->convert_to_alias_deriv(dF);
+	s->convert_to_alias_deriv(dF); // we need two factors of drho_dx
+
 	matrix_dot_v1(dF,ff);
+
+	s->square_and_scale_with_d2rho_dx2(ff_save);
+	dF += ff_save;
+
 	dF.set(ff);
 	s->endForceCalculation(); // Need to do this again to make sure that the boundaries are fixed, when demanded	
       }

diff --git a/src/FMT_Species.cpp b/src/FMT_Species.cpp
@@ -119,6 +119,31 @@ void FMT_Species::convert_to_alias_deriv(DFT_Vec &x, DFT_Vec &dF_dRho) const
   }
 }
 
+void FMT_Species::square_and_scale_with_d2rho_dx2(DFT_Vec &ff) const 
+{
+  long pos;
+  const double etamin = dmin*(4*M_PI*getHSD()*getHSD()*getHSD()/3);
+  const double c = (1.0-etamin)/density_->dV();
+
+  DFT_Vec x(ff.size());  
+  get_density_alias(x);  
+
+
+  #ifdef USE_OMP
+  #pragma omp parallel for  private(pos)  schedule(static)
+  #endif
+  for(pos=0;pos<x.size();pos++)
+  {
+    double y = x.get(pos);
+    double df = ff.get(pos);
+
+    double drho_dx   = 2*c*y/((1+y*y)*(1+y*y));
+    double d2rho_dx2 = (2*c/((1+y*y)*(1+y*y))) - 4*y*drho_dx/(1+y*y);
+
+    ff.set(pos, df*df*d2rho_dx2);
+  }
+}
+
 void FMT_Species::convert_to_alias_increment(DFT_Vec &dRho) const
 {
   DFT_Vec x; get_density_alias(x);