From 3355e82837b6a0476426fb8bc54d2428014fea58 Mon Sep 17 00:00:00 2001 From: Eric Bylaska Date: Thu, 20 Jul 2023 16:47:38 -0700 Subject: [PATCH] Fixed paw ortho bug when using processor groups. D1dB_SumAll was being used instead of D1dB_Vector_SumAll..EJB --- src/nwpw/nwpwlib/D1dB/D1dB-mpi.F | 14 +-- src/nwpw/nwpwlib/D1dB/Dne.F | 13 ++- src/nwpw/nwpwlib/control/control.F | 6 +- src/nwpw/nwpwlib/ion/ion.F | 12 +- src/nwpw/pspw/cpsd/GNUmakefile | 16 +-- src/nwpw/pspw/cpsd/inner_loop.F | 10 +- src/nwpw/pspw/cpsd/psi_lmbda_paw.f | 2 +- src/nwpw/pspw/cpsd/psi_lmbda_paw1.f | 164 ++++++++++++++++++++++++++++ 8 files changed, 200 insertions(+), 37 deletions(-) create mode 100644 src/nwpw/pspw/cpsd/psi_lmbda_paw1.f diff --git a/src/nwpw/nwpwlib/D1dB/D1dB-mpi.F b/src/nwpw/nwpwlib/D1dB/D1dB-mpi.F index 8f59f1f9a7..63639ed74c 100644 --- a/src/nwpw/nwpwlib/D1dB/D1dB-mpi.F +++ b/src/nwpw/nwpwlib/D1dB/D1dB-mpi.F @@ -4,11 +4,10 @@ * * *********************************** -* * * -* * D1dB_SumAll * -* * * +* * * +* * D1dB_SumAll * +* * * * *********************************** - subroutine D1dB_SumAll(sum) c implicit none real*8 sum @@ -191,11 +190,10 @@ c implicit none * *********************************** -* * * -* * D1dB_Vector_SumAll * -* * * +* * * +* * D1dB_Vector_SumAll * +* * * * *********************************** - subroutine D1dB_Vector_SumAll(n,sum) c implicit none integer n diff --git a/src/nwpw/nwpwlib/D1dB/Dne.F b/src/nwpw/nwpwlib/D1dB/Dne.F index d5a7bade9b..44c7617a93 100644 --- a/src/nwpw/nwpwlib/D1dB/Dne.F +++ b/src/nwpw/nwpwlib/D1dB/Dne.F @@ -3935,10 +3935,9 @@ c **************************************** c **************************************** c * * -c * Dneall_m_scale_s21 * +c * Dneall_m_scale_s21 * c * * c **************************************** - subroutine Dneall_m_scale_s21(mb,dte,s21) implicit none integer mb @@ -4008,9 +4007,10 @@ c **************************************** end + c **************************************** c * * -c * Dneall_m_scale_s11 * +c * Dneall_m_scale_s11 * c * * c **************************************** @@ -4507,7 +4507,7 @@ c **************************************** real*8 ddot external ddot - + if (parallelized) then if (mb.eq.0) then nn = ne(1) + ne(2) @@ -4545,8 +4545,8 @@ c **************************************** end do !$OMP END MASTER !$OMP BARRIER - call D1dB_SumAll(nn*nprj,dbl_mb(sw1b(1))) - call D1dB_SumAll(nn*nprj,dbl_mb(sw2b(1))) + call D1dB_Vector_SumAll(nn*nprj,dbl_mb(sw1b(1))) + call D1dB_Vector_SumAll(nn*nprj,dbl_mb(sw2b(1))) if (mparallelized) then call Parallel_taskid(taskid) @@ -4609,6 +4609,7 @@ c **************************************** if (.not.value) > call errquit('Dneall_m_add_sw1sw2:pop stack',1,MA_ERR) else + if (mb.eq.0) then ms1 = 1 ms2 = ispin diff --git a/src/nwpw/nwpwlib/control/control.F b/src/nwpw/nwpwlib/control/control.F index b790ba62df..869c2f45e5 100644 --- a/src/nwpw/nwpwlib/control/control.F +++ b/src/nwpw/nwpwlib/control/control.F @@ -3,9 +3,9 @@ * * *********************************** -* * * -* * control_read * -* * * +* * * +* * control_read * +* * * * *********************************** logical function control_read(code_in,rtdb) implicit none diff --git a/src/nwpw/nwpwlib/ion/ion.F b/src/nwpw/nwpwlib/ion/ion.F index f448546cae..a32dfa7d7f 100644 --- a/src/nwpw/nwpwlib/ion/ion.F +++ b/src/nwpw/nwpwlib/ion/ion.F @@ -1460,9 +1460,9 @@ c end if * *************************** -* * * -* * ion_chargeexist * -* * * +* * * +* * ion_chargeexist * +* * * * *************************** logical function ion_chargeexist() implicit none @@ -1475,9 +1475,9 @@ c end if end * *************************** -* * * -* * ion_nion * -* * * +* * * +* * ion_nion * +* * * * *************************** integer function ion_nion() implicit none diff --git a/src/nwpw/pspw/cpsd/GNUmakefile b/src/nwpw/pspw/cpsd/GNUmakefile index 1b2681be5b..02e02c44e3 100644 --- a/src/nwpw/pspw/cpsd/GNUmakefile +++ b/src/nwpw/pspw/cpsd/GNUmakefile @@ -4,13 +4,13 @@ OBJ_OPTIMIZE = \ cpsd_input.o \ cpmd_input.o \ - cpsdv5.o \ - cpmdv5.o \ + cpsdv5.o \ + cpmdv5.o \ mmsdv1.o \ - inner_loop_md.o \ - mm_loop.o \ - psi_lmbda.o psi_lmbda_sic.o psi_lmbda_paw.o psi_lmbda2.o psi_lmbda3.o \ - psi_lmbda_omp.o \ + inner_loop_md.o \ + mm_loop.o \ + psi_lmbda.o psi_lmbda_sic.o psi_lmbda_paw.o psi_lmbda2.o psi_lmbda3.o \ + psi_lmbda_omp.o psi_lmbda_paw1.o \ fei_output.o \ dipole_motion.o \ pspw_et.o \ @@ -23,9 +23,9 @@ LIBRARY = libpspw.a - USES_BLAS = psi_lmbda.f psi_lmbda2.f psi_lmbda3.f psi_lmbda_sic.f psi_lmbda_paw.f \ + USES_BLAS = psi_lmbda.f psi_lmbda2.f psi_lmbda3.f psi_lmbda_sic.f psi_lmbda_paw.f psi_lmbda_paw1.f \ inner_loop.F inner_loop_md.F cpsdv5.F cpmdv5.F psi_lmbda2.f psi_lmbda_omp.F \ - mmsdv1.F inner_loop_qmmm_step.F cpmd_qmmm_start.F cpmd_qmmm_stop.F mm_loop.F pspw_et.F + mmsdv1.F inner_loop_qmmm_step.F cpmd_qmmm_start.F cpmd_qmmm_stop.F mm_loop.F pspw_et.F LIB_DEFINES = diff --git a/src/nwpw/pspw/cpsd/inner_loop.F b/src/nwpw/pspw/cpsd/inner_loop.F index e6b46bedd2..0383a12b39 100644 --- a/src/nwpw/pspw/cpsd/inner_loop.F +++ b/src/nwpw/pspw/cpsd/inner_loop.F @@ -670,13 +670,13 @@ c call Pack_cc_Sum(1,psi2(1,n),psi1(1,n),psi2(1,n)) call Dneall_f_Sortho(ms,psi2,psi_r,npack1) end if else - call psi_lmbda_paw(ispin,ne,nemax,npack1,psi_r,psi2, - > dte, - > lmd,dbl_mb(tmp_L(1)),ierr) + call psi_lmbda_paw(ispin,ne,nemax,npack1,psi_r,psi2,dte, + > lmd,dbl_mb(tmp_L(1)),ierr) +c call psi_lmbda_paw1(ispin,ne,nemax,npack1,psi_r,psi2,dte, +c > lmd,dbl_mb(tmp_L(1)),ierr) end if else if (fractional) then - call psi_lmbda2(ispin,ne,nemax,npack1,psi1,psi2, - > dte,occ1, + call psi_lmbda2(ispin,ne,nemax,npack1,psi1,psi2,dte,occ1, > lmd,dbl_mb(tmp_L(1)),ierr) !else if (sic.or.hfx) then else if (sic) then diff --git a/src/nwpw/pspw/cpsd/psi_lmbda_paw.f b/src/nwpw/pspw/cpsd/psi_lmbda_paw.f index db73651967..7a2442cefb 100644 --- a/src/nwpw/pspw/cpsd/psi_lmbda_paw.f +++ b/src/nwpw/pspw/cpsd/psi_lmbda_paw.f @@ -2,7 +2,7 @@ c $Id: psi_lmbda_paw.f 21177 2011-10-10 17:09:43Z bylaska $ * *********************************** * * * -* * psi_lmbda_paw * +* * psi_lmbda_paw * * * * * *********************************** diff --git a/src/nwpw/pspw/cpsd/psi_lmbda_paw1.f b/src/nwpw/pspw/cpsd/psi_lmbda_paw1.f new file mode 100644 index 0000000000..8f8f748048 --- /dev/null +++ b/src/nwpw/pspw/cpsd/psi_lmbda_paw1.f @@ -0,0 +1,164 @@ +* +* $Id$ +* + +* ******************************************** +* * * +* * psi_lmbda_paw1 * +* * * +* ******************************************** + subroutine psi_lmbda_paw1(ispin,ne,nemaxq,npack1, + > psi1,psi2, + > dte, + > lmbda,tmp,ierr) + implicit none + integer ispin,ne(2),nemaxq,npack1 + complex*16 psi1(npack1,nemaxq) + complex*16 psi2(npack1,nemaxq) + real*8 dte + real*8 lmbda(*) + real*8 tmp(*) + integer ierr + + integer MASTER + parameter (MASTER=0) + +* **** local variables **** + logical notgram + integer taskid + integer ms,it + integer nn + integer s11,s12,s21,s22,st1,st2,sa1,sa0 + real*8 adiff + +* :::: iteration limit and tolerence for non-liner equations :::: + integer itrlmd + real*8 convg,convg2 + parameter (itrlmd=120, convg=1.0d-15,convg2=1.0e-12) + +* **** external functions **** + real*8 Dneall_m_dmax + external Dneall_m_dmax + integer Parallel_threadid + external Parallel_threadid + + call nwpw_timing_start(3) + +c call Dneall_m_size(1,nn) +c s11 = 0*nn + 1 +c s21 = 1*nn + 1 +c s22 = 2*nn + 1 +c s12 = 3*nn + 1 +c sa0 = 4*nn + 1 +c sa1 = 5*nn + 1 +c st1 = 6*nn + 1 +c st2 = 7*nn + 1 +c call dcopy(8*nn,0.0d0,0,tmp,1) + + +*:::::::::::::::::::::: Lagrangian multipliers :::::::::::::::::::::: + DO 640 ms=1,ispin + notgram=.true. + IF(ne(ms).le.0) GO TO 640 + + call Dneall_m_size(ms,nn) + s11 = 0*nn + 1 + s21 = 1*nn + 1 + s22 = 2*nn + 1 + s12 = 3*nn + 1 + sa0 = 4*nn + 1 + sa1 = 5*nn + 1 + st1 = 6*nn + 1 + st2 = 7*nn + 1 + +ccccccccc Uncomment in this block and then run ch3cl-pspw-sd.nw cccccccccccc + call Dneall_ffm_sym_Multiply(ms,psi2,psi2,npack1,tmp(s22)) + call Dneall_ffm_Multiply(ms,psi2,psi1,npack1,tmp(s21)) + call Dneall_ffm_Multiply(ms,psi1,psi2,npack1,tmp(s12)) + call Dneall_ffm_sym_Multiply(ms,psi1,psi1,npack1,tmp(s11)) + !call Dne_ffm_combo_sym_Multiply(ms,psi1,psi2,npack1,tmp(s11)) + +* *** add the extra paw overlap part *** + call psp_add_paw_extra_overlap1(ms,psi2,tmp(s22)) + call psp_add_paw_extra_overlap2(ms,psi2,psi1,tmp(s21)) + call psp_add_paw_extra_overlap2(ms,psi1,psi2,tmp(s12)) + call psp_add_paw_extra_overlap1(ms,psi1,tmp(s11)) + +ccccccccc Uncomment in this block and then run ch3cl-pspw-sd.nw cccccccccccc + +* ***** scale the overlap matrices **** + call Dneall_m_scale_s22(ms,dte,tmp(s22)) + call Dneall_m_scale_s21(ms,dte,tmp(s21)) + call Dneall_m_scale_s21(ms,dte,tmp(s12)) + call Dneall_m_scale_s11(ms,dte,tmp(s11)) + + !call dcopy(nn,tmp(s22),1,tmp(sa0),1) + call Parallel_shared_vector_copy(.true.,nn,tmp(s22),tmp(sa0)) + !call Parallel_shared_vector_copy(.true.,nn,tmp(s21),tmp(s12)) + + + do it=1,itrlmd + !call dcopy(nn,tmp(s22),1,tmp(sa1),1) + call Parallel_shared_vector_copy(.true.,nn,tmp(s22),tmp(sa1)) + + call Dneall_mmm_Multiply0(ms, + > tmp(s21),tmp(sa0),1.0d0, + > tmp(sa1),1.0d0) + call Dneall_mmm_Multiply0(ms, + > tmp(sa0),tmp(s12),1.0d0, + > tmp(sa1),1.0d0) + call Dneall_mmm_Multiply0(ms, + > tmp(s11),tmp(sa0),1.0d0, + > tmp(st1),0.0d0) + call Dneall_mmm_Multiply(ms, + > tmp(sa0),tmp(st1),1.0d0, + > tmp(sa1),1.0d0) + !call dcopy(nn,tmp(sa1),1,tmp(st1),1) + call Parallel_shared_vector_copy(.true.,nn,tmp(sa1),tmp(st1)) + call daxpy_omp(nn,(-1.0d0),tmp(sa0),1,tmp(st1),1) + + adiff = Dneall_m_dmax(ms,tmp(st1)) + + if(adiff.lt.convg) GO TO 630 + if (adiff.gt.1.0d10) go to 620 + !call dcopy(nn,tmp(sa1),1,tmp(sa0),1) + call Parallel_shared_vector_copy(.true.,nn,tmp(sa1),tmp(sa0)) + end do + if (adiff.lt.convg2) GO TO 630 + + 620 continue + ierr=10 + call Parallel_taskid(taskid) + if (taskid.eq.MASTER) then + WRITE(6,*) + > 'Warning: Lagrange Multiplier tolerance too high:',adiff + WRITE(6,*) ' +Try using a smaller time step' + WRITE(6,*) ' +Gram-Schmidt being performed, spin:',ms + end if +c call Dneall_f_ortho(ms,psi2,npack1) + call Dneall_f_GramSchmidt(ms,psi2,npack1) + notgram = .false. + + 630 continue + +* ::::::::::::::::: correction due to the constraint ::::::::::::::::: + if (notgram) then + call Dneall_fmf_Multiply(ms, + > psi1,npack1, + > tmp(sa1), dte, + > psi2,1.0d0) + end if + call Dneall_mm_Expand(ms,tmp(sa1),lmbda) + 640 continue + +c*::::::::::::::::: correction due to the constraint ::::::::::::::::: +c call Dneall_fmf_Multiply(0, +c > psi1,npack1, +c > lmbda, dte, +c > psi2,1.0d0) + + call nwpw_timing_end(3) + + return + end +