From 97d3c32dc348a40e0c56bc2dbd846603def2455f Mon Sep 17 00:00:00 2001 From: JFL Date: Thu, 7 Sep 2023 11:06:08 +0200 Subject: [PATCH] synch --- include/DFT_LinAlg.h | 1 + include/Density.h | 3 +- src/DFT.cpp | 109 +++++++++++++++++++++++-------------------- src/FMT.cpp | 6 ++- 4 files changed, 65 insertions(+), 54 deletions(-) diff --git a/include/DFT_LinAlg.h b/include/DFT_LinAlg.h index a65aae6..dd1971e 100644 --- a/include/DFT_LinAlg.h +++ b/include/DFT_LinAlg.h @@ -270,6 +270,7 @@ class DFT_FFT void MultBy(double val) { RealSpace_.MultBy(val); FourierSpace_.MultBy(val);} // doesn't change is_dirty_ void IncrementBy(DFT_Vec & v) { RealSpace_.IncrementBy(v); is_dirty_ = true;} + void DecrementBy(DFT_Vec & v) { RealSpace_.DecrementBy(v); is_dirty_ = true;} void IncrementBy(unsigned pos, double val) { RealSpace_.IncrementBy(pos,val); is_dirty_ = true;} friend ostream &operator<<(ostream &of, const DFT_FFT &v) diff --git a/include/Density.h b/include/Density.h index b348029..3d3cb40 100644 --- a/include/Density.h +++ b/include/Density.h @@ -100,7 +100,8 @@ class Density : public Lattice // do stuff to the density void operator*=(double a) { Density_.MultBy(a);} - void operator+=(DFT_Vec &v) { Density_.IncrementBy(v);} + void operator+=(DFT_Vec &v) { Density_.IncrementBy(v);} + void operator-=(DFT_Vec &v) { Density_.DecrementBy(v);} void shift(DFT_Vec &direction, double scale) { Density_.Real().IncrementBy_Scaled_Vector(direction, scale);} void center_cluster(bool homogeneousBoundary); diff --git a/src/DFT.cpp b/src/DFT.cpp index 40292ae..a837a7b 100644 --- a/src/DFT.cpp +++ b/src/DFT.cpp @@ -154,6 +154,12 @@ double DFT::calculateFreeEnergyAndDerivatives(bool onlyFex, bool H_dot_Force) rms_force_ += s->getDF().euclidean_norm(); } rms_force_ /= sqrt(allSpecies_.size()*get_Ntot()); + + // The following is used when we want to minimize L = (dF/dx_I)^2 = (F_I rho(x_I)')^2 = (F_I)^2 (rho(x_I)')^2 + // In this case, the force is + // 0.5*dL/dx_J = F_I (rho(x_I)')^2 F_IJ rho(x_J)' + (F_J)^2 (2*rho(x_J)'*rho(x_J)'') + // = [F_I (rho(x_I)')^2 F_IJ + (F_I)^2 (2*rho(x_J)'')]*rho(x_J)' + // Here, we calculate the term in square brackets because the last factor is added in the minimization routines. if(H_dot_Force) for(auto &s: allSpecies_) @@ -161,19 +167,19 @@ double DFT::calculateFreeEnergyAndDerivatives(bool onlyFex, bool H_dot_Force) if(s->is_mass_fixed()) throw std::runtime_error("Makes no sense to do H_dot_Force when mass is fixed"); DFT_Vec &dF = s->getDF(); - DFT_Vec ff(dF.size()); - - DFT_Vec ff_save(dF); + DFT_Vec dF_copy(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->convert_to_alias_deriv(dF); // we need two factors of drho_dx + + DFT_Vec first_term(dF.size()); + matrix_dot_v1(dF,first_term); - s->square_and_scale_with_d2rho_dx2(ff_save); - dF += ff_save; + s->square_and_scale_with_d2rho_dx2(dF_copy); + dF_copy *= 2; + first_term += dF_copy; - dF.set(ff); + dF.set(first_term); s->endForceCalculation(); // Need to do this again to make sure that the boundaries are fixed, when demanded } @@ -285,14 +291,14 @@ void DFT::matrix_dot_v_intern(const vector &v_in, vector &resu vector v(allSpecies_.size()); for(int s=0;sconvert_to_density_increment(vv); + { + DFT_Vec vv = v_in[s].cReal(); + if (is_matrix_in_alias_coordinates()) allSpecies_[s]->convert_to_density_increment(vv); - v[s].initialize(Nx, Ny, Nz); - v[s].Real().set(vv); - v[s].do_real_2_fourier(); - } + v[s].initialize(Nx, Ny, Nz); + v[s].Real().set(vv); + v[s].do_real_2_fourier(); + } // I would like to do this but it violates the const declaration of v // Well this is now implemented in the lines above ... @@ -303,58 +309,59 @@ void DFT::matrix_dot_v_intern(const vector &v_in, vector &resu for(int s=0;sis_fixed_boundary()) for(unsigned p=0;pgetLattice().get_Nboundary();p++) - { - unsigned pp = allSpecies_[s]->getLattice().boundary_pos_2_pos(p); - if(fabs(v[s].cReal().get(pp)) > 0.0) - throw std::runtime_error("Input vector v must have zero boundary entries in DFT::hessian_dot_v when the species has fixed boundaries"); - } + { + unsigned pp = allSpecies_[s]->getLattice().boundary_pos_2_pos(p); + if(fabs(v[s].cReal().get(pp)) > 0.0) + throw std::runtime_error("Input vector v must have zero boundary entries in DFT::hessian_dot_v when the species has fixed boundaries"); + } // ideal gas contribution: v_i/n_i if(full_hessian_) - { - double dV = allSpecies_[0]->getDensity().dV(); + { + double dV = allSpecies_[0]->getDensity().dV(); - for(int s=0;sgetDensity().get(pos)); - } + for(int s=0;sgetDensity().get(pos)); + } // Hard-sphere if(fmt_) - { - try {fmt_->add_second_derivative(v,result, allSpecies_);} - catch( Eta_Too_Large_Exception &e) {throw e;} - } else for(auto &species : allSpecies_) // needed for interactions - its done automatically if there is an fmt evaluation - species->doFFT(); + { + try {fmt_->add_second_derivative(v,result, allSpecies_);} + catch( Eta_Too_Large_Exception &e) {throw e;} + } else for(auto &species : allSpecies_) // needed for interactions - its done automatically if there is an fmt evaluation + species->doFFT(); // Mean field for(auto &interaction: DFT::Interactions_) interaction->add_second_derivative(v,result); - if (is_matrix_in_alias_coordinates()) for(int s=0;sgetDF(); - DFT_Vec d2rho; allSpecies_[s]->get_second_derivatives_of_density_wrt_alias(d2rho); - df_term.Schur(df_term, d2rho); - df_term.Schur(df_term, v_in[s].cReal()); // Note that we use here the vector in alias coordinates + if (is_matrix_in_alias_coordinates()) + for(int s=0;sgetDF(); + DFT_Vec d2rho; allSpecies_[s]->get_second_derivatives_of_density_wrt_alias(d2rho); + df_term.Schur(df_term, d2rho); + df_term.Schur(df_term, v_in[s].cReal()); // Note that we use here the vector in alias coordinates - // Add all terms together - allSpecies_[s]->convert_to_alias_deriv(result[s]); - result[s].IncrementBy(df_term); - } + // Add all terms together + allSpecies_[s]->convert_to_alias_deriv(result[s]); + result[s].IncrementBy(df_term); + } // Remove boundary terms if the boundary is fixed for(int s=0;sis_fixed_boundary()) - { - long p = 0; - do{result[s].set(p,0.0);} while(get_next_boundary_point(p)); - } + { + long p = 0; + do{result[s].set(p,0.0);} while(get_next_boundary_point(p)); + } } // computes result_I = F_{I I+J} for fixed J diff --git a/src/FMT.cpp b/src/FMT.cpp index 699a046..46f898b 100644 --- a/src/FMT.cpp +++ b/src/FMT.cpp @@ -592,7 +592,8 @@ void FMT::add_second_derivative(int jx, int jy, int jz, const vector & for(int a = 0;a & result.do_fourier_2_real(); d2F[0].IncrementBy_Scaled_Vector(result.Real(),dV/Ntot); // Ntot because FFTW normalization - } + } + cout << endl; }