@@ -86,28 +86,18 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
8686 if (isforce)
8787 {
8888 fcs.create (nat, 3 );
89- foverlap.create (nat, 3 );
90- ftvnl_dphi.create (nat, 3 );
91- fvnl_dbeta.create (nat, 3 );
92- fvl_dvl.create (nat, 3 );
93- fewalds.create (nat, 3 );
94- fcc.create (nat, 3 );
95- fscc.create (nat, 3 );
89+ foverlap.create (nat, 3 ); // overlap force
90+ ftvnl_dphi.create (nat, 3 ); // pulay force of NAO
91+ fvnl_dbeta.create (nat, 3 ); // pulay force of non-local projectors
92+ fvl_dvl.create (nat, 3 ); // force from local potentials
93+ fewalds.create (nat, 3 ); // Ewald force
94+ fcc.create (nat, 3 ); // force due to core correction
95+ fscc.create (nat, 3 ); // force due to self-consistent field
9696 fvnl_dalpha.create (nat, 3 ); // deepks
9797
9898 // calculate basic terms in Force, same method with PW base
99- this ->calForcePwPart (ucell,
100- fvl_dvl,
101- fewalds,
102- fcc,
103- fscc,
104- pelec->f_en .etxc ,
105- pelec->vnew ,
106- pelec->vnew_exist ,
107- pelec->charge ,
108- rhopw,
109- locpp,
110- sf);
99+ this ->calForcePwPart (ucell, fvl_dvl, fewalds, fcc, fscc, pelec->f_en .etxc ,
100+ pelec->vnew , pelec->vnew_exist , pelec->charge , rhopw, locpp, sf);
111101 }
112102
113103 // total stress : ModuleBase::matrix scs
@@ -139,63 +129,30 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
139129 svnl_dalpha.create (3 , 3 );
140130
141131 // calculate basic terms in Stress, similar method with PW base
142- this ->calStressPwPart (ucell,
143- sigmadvl,
144- sigmahar,
145- sigmaewa,
146- sigmacc,
147- sigmaxc,
148- pelec->f_en .etxc ,
149- pelec->charge ,
150- rhopw,
151- locpp,
152- sf);
132+ this ->calStressPwPart (ucell, sigmadvl, sigmahar, sigmaewa, sigmacc,
133+ sigmaxc, pelec->f_en .etxc , pelec->charge , rhopw, locpp, sf);
153134 }
154135
155136 // ! atomic forces from integration (4 terms)
156- this ->integral_part (PARAM.globalv .gamma_only_local ,
157- isforce,
158- isstress,
159- ucell,
160- gd,
161- fsr,
162- pelec,
163- psi,
164- foverlap,
165- ftvnl_dphi,
166- fvnl_dbeta,
167- fvl_dphi,
168- soverlap,
169- stvnl_dphi,
170- svnl_dbeta,
171- svl_dphi,
172- #ifdef __MLALGO
173- fvnl_dalpha,
174- svnl_dalpha,
175- deepks.ld ,
176- #endif
177- gint_gamma,
178- gint_k,
179- two_center_bundle,
180- orb,
181- pv,
182- kv);
137+ this ->integral_part (PARAM.globalv .gamma_only_local , isforce, isstress,
138+ ucell, gd, fsr, pelec, psi, foverlap, ftvnl_dphi,
139+ fvnl_dbeta, fvl_dphi, soverlap, stvnl_dphi, svnl_dbeta,
140+ svl_dphi, fvnl_dalpha, svnl_dalpha, deepks, gint_gamma,
141+ gint_k, two_center_bundle, orb, pv, kv);
142+
183143 // calculate force and stress for Nonlocal part
184144 if (PARAM.inp .nspin == 1 || PARAM.inp .nspin == 2 )
185145 {
186146 hamilt::NonlocalNew<hamilt::OperatorLCAO<T, double >> tmp_nonlocal (nullptr ,
187- kv.kvec_d ,
188- nullptr ,
189- &ucell,
190- orb.cutoffs (),
191- &gd,
192- two_center_bundle.overlap_orb_beta .get ());
147+ kv.kvec_d , nullptr , &ucell, orb.cutoffs (), &gd, two_center_bundle.overlap_orb_beta .get ());
193148
194149 const auto * dm_p = dynamic_cast <const elecstate::ElecStateLCAO<T>*>(pelec)->get_DM ();
150+
195151 if (PARAM.inp .nspin == 2 )
196152 {
197153 const_cast <elecstate::DensityMatrix<T, double >*>(dm_p)->switch_dmr (1 );
198154 }
155+
199156 const hamilt::HContainer<double >* dmr = dm_p->get_DMR_pointer (1 );
200157 tmp_nonlocal.cal_force_stress (isforce, isstress, dmr, fvnl_dbeta, svnl_dbeta);
201158 if (PARAM.inp .nspin == 2 )
@@ -206,12 +163,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
206163 else if (PARAM.inp .nspin == 4 )
207164 {
208165 hamilt::NonlocalNew<hamilt::OperatorLCAO<std::complex <double >, std::complex <double >>> tmp_nonlocal (
209- nullptr ,
210- kv.kvec_d ,
211- nullptr ,
212- &ucell,
213- orb.cutoffs (),
214- &gd,
166+ nullptr , kv.kvec_d , nullptr , &ucell, orb.cutoffs (), &gd,
215167 two_center_bundle.overlap_orb_beta .get ());
216168
217169 // calculate temporary complex DMR for nonlocal force&stress
@@ -777,11 +729,9 @@ void Force_Stress_LCAO<double>::integral_part(const bool isGammaOnly,
777729 ModuleBase::matrix& stvnl_dphi,
778730 ModuleBase::matrix& svnl_dbeta,
779731 ModuleBase::matrix& svl_dphi,
780- #if __MLALGO
781732 ModuleBase::matrix& fvnl_dalpha,
782733 ModuleBase::matrix& svnl_dalpha,
783- LCAO_Deepks<double >& ld,
784- #endif
734+ Setup_DeePKS<double >& deepks,
785735 Gint_Gamma& gint_gamma, // mohan add 2024-04-01
786736 Gint_k& gint_k, // mohan add 2024-04-01
787737 const TwoCenterBundle& two_center_bundle,
@@ -790,30 +740,11 @@ void Force_Stress_LCAO<double>::integral_part(const bool isGammaOnly,
790740 const K_Vectors& kv)
791741{
792742
793- flk.ftable (isforce,
794- isstress,
795- fsr, // mohan add 2024-06-15
796- ucell,
797- gd,
798- psi,
799- pelec,
800- foverlap,
801- ftvnl_dphi,
802- fvnl_dbeta,
803- fvl_dphi,
804- soverlap,
805- stvnl_dphi,
806- svnl_dbeta,
807- svl_dphi,
808- #if __MLALGO
809- fvnl_dalpha,
810- svnl_dalpha,
811- ld,
812- #endif
813- gint_gamma,
814- two_center_bundle,
815- orb,
816- pv);
743+ flk.ftable (isforce, isstress, fsr, ucell, gd, psi, pelec,
744+ foverlap, ftvnl_dphi, fvnl_dbeta, fvl_dphi,
745+ soverlap, stvnl_dphi, svnl_dbeta, svl_dphi,
746+ fvnl_dalpha, svnl_dalpha, deepks, gint_gamma,
747+ two_center_bundle, orb, pv);
817748 return ;
818749}
819750
@@ -834,44 +765,21 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(const bool isGammaOn
834765 ModuleBase::matrix& stvnl_dphi,
835766 ModuleBase::matrix& svnl_dbeta,
836767 ModuleBase::matrix& svl_dphi,
837- #if __MLALGO
838768 ModuleBase::matrix& fvnl_dalpha,
839769 ModuleBase::matrix& svnl_dalpha,
840- LCAO_Deepks<std::complex <double >>& ld,
841- #endif
770+ Setup_DeePKS<std::complex <double >>& deepks,
842771 Gint_Gamma& gint_gamma,
843772 Gint_k& gint_k,
844773 const TwoCenterBundle& two_center_bundle,
845774 const LCAO_Orbitals& orb,
846775 const Parallel_Orbitals& pv,
847776 const K_Vectors& kv)
848777{
849- flk.ftable (isforce,
850- isstress,
851- fsr, // mohan add 2024-06-16
852- ucell,
853- gd,
854- psi,
855- pelec,
856- foverlap,
857- ftvnl_dphi,
858- fvnl_dbeta,
859- fvl_dphi,
860- soverlap,
861- stvnl_dphi,
862- svnl_dbeta,
863- svl_dphi,
864- #if __MLALGO
865- fvnl_dalpha,
866- svnl_dalpha,
867- ld,
868- #endif
869- gint_k,
870- two_center_bundle,
871- orb,
872- pv,
873- &kv,
874- this ->RA );
778+ flk.ftable (isforce, isstress, fsr, ucell, gd, psi, pelec,
779+ foverlap, ftvnl_dphi, fvnl_dbeta, fvl_dphi,
780+ soverlap, stvnl_dphi, svnl_dbeta, svl_dphi,
781+ fvnl_dalpha, svnl_dalpha, deepks, gint_k,
782+ two_center_bundle, orb, pv, &kv, this ->RA );
875783 return ;
876784}
877785
@@ -890,30 +798,20 @@ void Force_Stress_LCAO<T>::calStressPwPart(UnitCell& ucell,
890798 const Structure_Factor& sf)
891799{
892800 ModuleBase::TITLE (" Force_Stress_LCAO" , " calStressPwPart" );
893- // --------------------------------------------------------
801+
894802 // local pseudopotential stress:
895- // use charge density; plane wave; local pseudopotential;
896- // --------------------------------------------------------
897803 sc_pw.stress_loc (ucell, sigmadvl, rhopw, locpp.vloc , &sf, 0 , chr);
898804
899- // --------------------------------------------------------
900805 // hartree term
901- // --------------------------------------------------------
902806 sc_pw.stress_har (ucell, sigmahar, rhopw, 0 , chr);
903807
904- // --------------------------------------------------------
905808 // ewald stress: use plane wave only.
906- // --------------------------------------------------------
907809 sc_pw.stress_ewa (ucell, sigmaewa, rhopw, 0 ); // remain problem
908810
909- // --------------------------------------------------------
910811 // stress due to core correlation.
911- // --------------------------------------------------------
912812 sc_pw.stress_cc (sigmacc, rhopw, ucell, &sf, 0 , locpp.numeric , chr);
913813
914- // --------------------------------------------------------
915814 // stress due to self-consistent charge.
916- // --------------------------------------------------------
917815 for (int i = 0 ; i < 3 ; i++)
918816 {
919817 sigmaxc (i, i) = -etxc / ucell.omega ;
@@ -932,21 +830,9 @@ void Force_Stress_LCAO<T>::forceSymmetry(const UnitCell& ucell, ModuleBase::matr
932830 double d1, d2, d3;
933831 for (int iat = 0 ; iat < ucell.nat ; iat++)
934832 {
935- ModuleBase::Mathzone::Cartesian_to_Direct (fcs (iat, 0 ),
936- fcs (iat, 1 ),
937- fcs (iat, 2 ),
938- ucell.a1 .x ,
939- ucell.a1 .y ,
940- ucell.a1 .z ,
941- ucell.a2 .x ,
942- ucell.a2 .y ,
943- ucell.a2 .z ,
944- ucell.a3 .x ,
945- ucell.a3 .y ,
946- ucell.a3 .z ,
947- d1,
948- d2,
949- d3);
833+ ModuleBase::Mathzone::Cartesian_to_Direct (fcs (iat, 0 ), fcs (iat, 1 ), fcs (iat, 2 ),
834+ ucell.a1 .x , ucell.a1 .y , ucell.a1 .z , ucell.a2 .x , ucell.a2 .y , ucell.a2 .z ,
835+ ucell.a3 .x , ucell.a3 .y , ucell.a3 .z , d1, d2, d3);
950836
951837 fcs (iat, 0 ) = d1;
952838 fcs (iat, 1 ) = d2;
@@ -955,21 +841,9 @@ void Force_Stress_LCAO<T>::forceSymmetry(const UnitCell& ucell, ModuleBase::matr
955841 symm->symmetrize_vec3_nat (fcs.c );
956842 for (int iat = 0 ; iat < ucell.nat ; iat++)
957843 {
958- ModuleBase::Mathzone::Direct_to_Cartesian (fcs (iat, 0 ),
959- fcs (iat, 1 ),
960- fcs (iat, 2 ),
961- ucell.a1 .x ,
962- ucell.a1 .y ,
963- ucell.a1 .z ,
964- ucell.a2 .x ,
965- ucell.a2 .y ,
966- ucell.a2 .z ,
967- ucell.a3 .x ,
968- ucell.a3 .y ,
969- ucell.a3 .z ,
970- d1,
971- d2,
972- d3);
844+ ModuleBase::Mathzone::Direct_to_Cartesian (fcs (iat, 0 ), fcs (iat, 1 ), fcs (iat, 2 ),
845+ ucell.a1 .x , ucell.a1 .y , ucell.a1 .z , ucell.a2 .x , ucell.a2 .y , ucell.a2 .z ,
846+ ucell.a3 .x , ucell.a3 .y , ucell.a3 .z , d1, d2, d3);
973847
974848 fcs (iat, 0 ) = d1;
975849 fcs (iat, 1 ) = d2;
0 commit comments