Skip to content

change vdwd3 from Fortran to C++ #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 2 commits into from
Closed

change vdwd3 from Fortran to C++ #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d1fe94c
1.change vdwd3 fortran part to C++(parameter `vdw_abc` temporarily in…
1041176461 Apr 26, 2021
4f0f3f2
add three body part of vdwd3, parameter `vdw_abc` can be used now
1041176461 Apr 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions ABACUS.develop/source/Makefile.Objects
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@ xc_2.o\
xc_3.o \
vdwd2.o\
vdwd2_parameters.o\
dftd3_subroutine.o\
vdwd3.o\
vdwd3_parameters.o\
pseudopot_upf.o\
Expand Down Expand Up @@ -350,4 +349,4 @@ lsame.o\
psort_w.o\
s_cmp.o\
s_copy.o\
xerbla.o\
xerbla.o\
2 changes: 1 addition & 1 deletion ABACUS.develop/source/Makefile.system
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ CEREAL_INCLUDE_DIR = ${CEREAL_DIR}/include
# comment out by mohan 2021-02-06
#LIBS = -lifcore -lm -lpthread ${LAPACK_LIB} ${FFTW_LIB} ${ELPA_LIB} ${LIBXC_LIB}
LIBS = -lifcore -lm -lpthread ${LAPACK_LIB} ${FFTW_LIB} ${ELPA_LIB}
#LIBS = -liomp5 -lpthread -lm -ldl ${BOOST_LIB} ${LAPACK_LIB} ${FFTW_LIB} ${LPA_LIB} ${LIBXC_LIB}
#LIBS = -liomp5 -lpthread -lm -ldl ${BOOST_LIB} ${LAPACK_LIB} ${FFTW_LIB} ${ELPA_LIB} ${LIBXC_LIB}

INCLUDES = -I. -Icommands -I${BOOST_INCLUDE_DIR} -I${LAPACK_INCLUDE_DIR} -I${FFTW_INCLUDE_DIR} -I${LIBXC_INCLUDE_DIR} -I${CEREAL_INCLUDE_DIR}

Expand Down
60 changes: 19 additions & 41 deletions ABACUS.develop/source/input_conv.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
#include "src_pw/efield.h"
#include "src_lcao/global_fp.h"
#include "src_pw/vdwd2_parameters.h"
#include "src_pw/vdwd3.h"
#include "src_pw/vdwd3_parameters.h"
#include "src_io/chi0_hilbert.h"
#include "src_io/chi0_standard.h"
#include "src_io/epsilon0_pwscf.h"
Expand Down Expand Up @@ -213,46 +213,24 @@ void Input_Conv::Convert(void)
}
if(INPUT.vdw_method=="d3_0" || INPUT.vdw_method=="d3_bj")
{
vdwd3.vdwD3 = true;
vdwd3.s6 = std::stod(INPUT.vdw_s6);
vdwd3.s18 = std::stod(INPUT.vdw_s8);
vdwd3.rs6 = std::stod(INPUT.vdw_a1);
vdwd3.rs18 = std::stod(INPUT.vdw_a2);
vdwd3.abc = INPUT.vdw_abc;
if(INPUT.vdw_method=="d3_0")
{
vdwd3.version = 1;
}
if(INPUT.vdw_method=="d3_bj")
{
vdwd3.version = 2;
}
vdwd3.model = INPUT.vdw_model;
if(INPUT.vdw_model=="radius")
{
if(INPUT.vdw_radius_unit=="Bohr")
{
vdwd3.rthr2 = pow(std::stod(INPUT.vdw_radius),2);
}
else
{
vdwd3.rthr2 = pow((std::stod(INPUT.vdw_radius) * BOHR_TO_A),2);
}
if(INPUT.vdw_cn_thr_unit=="Bohr")
{
vdwd3.cn_thr2 = pow(INPUT.vdw_cn_thr,2);
}
else
{
vdwd3.cn_thr2 = pow((INPUT.vdw_cn_thr * BOHR_TO_A),2);
}
}
else if(INPUT.vdw_model=="period")
{
vdwd3.rep_vdw[0] = INPUT.vdw_period.x;
vdwd3.rep_vdw[1] = INPUT.vdw_period.y;
vdwd3.rep_vdw[2] = INPUT.vdw_period.z;
}
vdwd3_para.flag_vdwd3 = true;
vdwd3_para.s6 = std::stod(INPUT.vdw_s6);
vdwd3_para.s18 = std::stod(INPUT.vdw_s8);
vdwd3_para.rs6 = std::stod(INPUT.vdw_a1);
vdwd3_para.rs18 = std::stod(INPUT.vdw_a2);
vdwd3_para.abc = INPUT.vdw_abc;
vdwd3_para.version = INPUT.vdw_method;
vdwd3_para.model = INPUT.vdw_model;
if(INPUT.vdw_radius_unit=="Bohr")
vdwd3_para.rthr2 = pow(std::stod(INPUT.vdw_radius),2);
else
vdwd3_para.rthr2 = pow((std::stod(INPUT.vdw_radius) * BOHR_TO_A),2);
if(INPUT.vdw_cn_thr_unit=="Bohr")
vdwd3_para.cn_thr2 = pow(INPUT.vdw_cn_thr,2);
else
vdwd3_para.cn_thr2 = pow((INPUT.vdw_cn_thr * BOHR_TO_A),2);
if(INPUT.vdw_model=="period")
vdwd3_para.period = INPUT.vdw_period;
}

//----------------------------------------------------------
Expand Down
2 changes: 1 addition & 1 deletion ABACUS.develop/source/src_global/matrix3.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ Matrix3 operator -(const Matrix3 &m1, const Matrix3 &m2); //m1-m2
Matrix3 operator /(const Matrix3 &m,const double &s); //m/s
Matrix3 operator *(const Matrix3 &m1,const Matrix3 &m2); //m1*m2
Matrix3 operator *(const Matrix3 &m,const double &s); //m*s
Matrix3 operator *(double &s, const Matrix3 &m); //s*m
Matrix3 operator *(const double &s, const Matrix3 &m); //s*m
template<typename T> Vector3<double> operator *(const Matrix3 &m, const Vector3<T> &u); //m*u // Peize Lin change Vector3<T> 2017-01-10
template<typename T> Vector3<double> operator *(const Vector3<T> &u, const Matrix3 &m); //u*m // Peize Lin change Vector3<T> 2017-01-10

Expand Down
36 changes: 26 additions & 10 deletions ABACUS.develop/source/src_lcao/FORCE_STRESS.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
// new
#include "../src_pw/H_XC_pw.h"
#include "src_pw/vdwd2.h"
#include "src_pw/vdwd3.h"

double Force_Stress_LCAO::force_invalid_threshold_ev = 0.00;
double Force_Stress_LCAO::output_acc = 1.0e-8;
Expand Down Expand Up @@ -154,13 +155,28 @@ void Force_Stress_LCAO::getForceStress(
stress_vdw = vdwd2.get_stress().to_matrix();
}
}
// jiyy add 2019-05-18
else if(vdwd3.vdwD3)
// jiyy add 2019-05-18, update 2021-04-25
else if(vdwd3_para.flag_vdwd3)
{
if(isforce) force_vdw.create(nat,3);
if(isstress) stress_vdw.create(3,3);
vdwd3.force(isforce, isstress, force_vdw, stress_vdw);
}
if(isforce)
{
force_vdw.create(nat,3);
Vdwd3 vdwd3(ucell,vdwd3_para);
vdwd3.cal_force();
for(int iat=0; iat<ucell.nat; ++iat)
{
force_vdw(iat,0) = vdwd3.get_force()[iat].x;
force_vdw(iat,1) = vdwd3.get_force()[iat].y;
force_vdw(iat,2) = vdwd3.get_force()[iat].z;
}
}
if(isstress)
{
Vdwd3 vdwd3(ucell,vdwd3_para);
vdwd3.cal_stress();
stress_vdw = vdwd3.get_stress().to_matrix();
}
}
//implement force from E-field
matrix fefield;
if(EFIELD&&isforce)
Expand Down Expand Up @@ -228,7 +244,7 @@ void Force_Stress_LCAO::getForceStress(
fcs(iat, i) += force_dftu(iat, i);
}
//VDW force of vdwd2 or vdwd3
if(vdwd2_para.flag_vdwd2||vdwd3.vdwD3)
if(vdwd2_para.flag_vdwd2||vdwd3_para.flag_vdwd3)
{
fcs(iat,i) += force_vdw(iat,i);
}
Expand Down Expand Up @@ -306,7 +322,7 @@ void Force_Stress_LCAO::getForceStress(
{
this->print_force("EFIELD FORCE",fefield,1,ry);
}
if(vdwd2_para.flag_vdwd2||vdwd3.vdwD3)
if(vdwd2_para.flag_vdwd2||vdwd3_para.flag_vdwd3)
{
this->print_force("VDW FORCE",force_vdw,1,ry);
}
Expand Down Expand Up @@ -358,7 +374,7 @@ void Force_Stress_LCAO::getForceStress(
+ sigmahar(i,j);// hartree stress

//VDW stress from linpz and jiyy
if(vdwd2_para.flag_vdwd2||vdwd3.vdwD3)
if(vdwd2_para.flag_vdwd2||vdwd3_para.flag_vdwd3)
{
scs(i,j) += stress_vdw(i , j);
}
Expand Down Expand Up @@ -414,7 +430,7 @@ void Force_Stress_LCAO::getForceStress(
sc_pw.print_stress("cc STRESS",sigmacc,TEST_STRESS,ry);
// sc_pw.print_stress("NLCC STRESS",sigmacc,TEST_STRESS,ry);
sc_pw.print_stress("XC STRESS",sigmaxc,TEST_STRESS,ry);
if(vdwd2_para.flag_vdwd2||vdwd3.vdwD3)
if(vdwd2_para.flag_vdwd2||vdwd3_para.flag_vdwd3)
{
sc_pw.print_stress("VDW STRESS",sigmaxc,TEST_STRESS,ry);
}
Expand Down
10 changes: 6 additions & 4 deletions ABACUS.develop/source/src_lcao/LOOP_elec.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "src_ri/exx_abfs.h"
#include "src_ri/exx_opt_orb.h"
#include "src_pw/vdwd2.h"
#include "src_pw/vdwd3.h"


void LOOP_elec::solve_elec_stru(const int &istep)
Expand Down Expand Up @@ -147,11 +148,12 @@ void LOOP_elec::solver(const int &istep)
vdwd2.cal_energy();
en.evdw = vdwd2.get_energy();
}
// jiyy add 2019-05-18
else if(vdwd3.vdwD3)
// jiyy add 2019-05-18, 2021-04-25
else if(vdwd3_para.flag_vdwd3)
{
vdwd3.energy();
en.evdw = vdwd3.energy_result;
Vdwd3 vdwd3(ucell,vdwd3_para);
vdwd3.cal_energy();
en.evdw = vdwd3.get_energy();
}

// self consistent calculations for electronic ground state
Expand Down
8 changes: 5 additions & 3 deletions ABACUS.develop/source/src_lcao/LOOP_ions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include "ELEC_scf.h"
#include "src_global/sltk_atom_arrange.h"
#include "src_pw/vdwd2.h"
#include "src_pw/vdwd3.h"
#include "LCAO_descriptor.h"

LOOP_ions::LOOP_ions()
Expand Down Expand Up @@ -447,10 +448,11 @@ void LOOP_ions::final_scf(void)
vdwd2.cal_energy();
en.evdw = vdwd2.get_energy();
}
else if(vdwd3.vdwD3) //jiyy add 2019-05-18
else if(vdwd3_para.flag_vdwd3) //jiyy add 2019-05-18, update 2021-04-25
{
vdwd3.energy();
en.evdw = vdwd3.energy_result;
Vdwd3 vdwd3(ucell,vdwd3_para);
vdwd3.cal_energy();
en.evdw = vdwd3.get_energy();;
}


Expand Down
8 changes: 5 additions & 3 deletions ABACUS.develop/source/src_lcao/run_md.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include "ELEC_scf.h"
#include "src_global/sltk_atom_arrange.h"
#include "src_pw/vdwd2.h"
#include "src_pw/vdwd3.h"

Run_MD::Run_MD()
{}
Expand Down Expand Up @@ -447,10 +448,11 @@ void Run_MD::final_scf(void)
vdwd2.cal_energy();
en.evdw = vdwd2.get_energy();
}
else if(vdwd3.vdwD3) //jiyy add 2019-05-18
else if(vdwd3_para.flag_vdwd3) //jiyy add 2019-05-18, update 2021-04-25
{
vdwd3.energy();
en.evdw = vdwd3.energy_result;
Vdwd3 vdwd3(ucell,vdwd3_para);
vdwd3.cal_energy();
en.evdw = vdwd3.get_energy();
}

ELEC_scf es;
Expand Down
Loading