diff --git a/src/nwpw/band/minimizer/GNUmakefile b/src/nwpw/band/minimizer/GNUmakefile index 60fcbf8a29..a1a209897f 100644 --- a/src/nwpw/band/minimizer/GNUmakefile +++ b/src/nwpw/band/minimizer/GNUmakefile @@ -12,6 +12,7 @@ c_cgminimize2.o \ c_cgminimizeABC.o \ c_bfgsminimize.o \ + c_bfgsminimizeABC.o \ c_bybminimize0.o \ c_bybminimize.o \ c_bybminimize2.o \ diff --git a/src/nwpw/band/minimizer/c_bfgsminimizeABC.F b/src/nwpw/band/minimizer/c_bfgsminimizeABC.F new file mode 100644 index 0000000000..7d28aaba16 --- /dev/null +++ b/src/nwpw/band/minimizer/c_bfgsminimizeABC.F @@ -0,0 +1,280 @@ +c +c $Id$ +c + +* ************************************************************ +* * * +* * Limited Memory BFGS ABC routine * +* * * +* ************************************************************ + + subroutine c_bfgsminimizeABC(E,deltae,deltac, + > current_iteration, + > minimizer) + implicit none + real*8 E(*) + real*8 deltae,deltac + integer current_iteration + integer minimizer + +#include "bafdecls.fh" +#include "errquit.fh" + +* **** local variables **** + + real*8 deltat_min + parameter (deltat_min=1.0d-2) + + integer G0_tag,S0_tag,G0_shift,S0_shift,G1_tag,G2_tag + real*8 E0,dE0,Ermdr + + logical precondition + common / cgsd_block2 / precondition + + real*8 tole,tolc + real*8 ehartree,eorbit,exc,pxc,eion + real*8 Enew,Eold,Estart + common / c_cgsd_block / Enew,Eold,Estart + + real*8 tmin,dte,sum0,sum1 + common / c_bfgs_block / tmin,dte,sum0,sum1 + + integer it,it_in + real*8 deltat + real*8 max_sigma,dt + real*8 ehfx,phfx + + integer ispin,neq(2),neall,npack1,nbrillq + +* **** define pointers **** + integer Y,U,HY,HU + +* **** external functions **** + integer control_it_in,cpsi_neq,cpsi_ispin,cpsi_nbrillq + integer cpsi_data_push_stack,cpsi_data_get_allptr + real*8 control_tole,control_tolc + real*8 cpsi_geodesic_energy + real*8 cpsi_geodesic_denergy + real*8 c_rho_error + real*8 c_dng_1ehartree + real*8 cpsi_1ke + real*8 cpsi_1vl + real*8 cpsi_1vnl + real*8 cpsi_1vnlso + real*8 c_rho_1exc + real*8 c_rho_1pxc + real*8 ewald_e + real*8 cpsi_1eorbit + real*8 linesearch + + external control_it_in,cpsi_neq,cpsi_ispin,cpsi_nbrillq + external cpsi_data_push_stack,cpsi_data_get_allptr + external control_tole,control_tolc + external cpsi_geodesic_energy + external cpsi_geodesic_denergy + external c_rho_error + external c_dng_1ehartree + external cpsi_1ke + external cpsi_1vl + external cpsi_1vnl + external cpsi_1vnlso + external c_rho_1exc + external c_rho_1pxc + external ewald_e + external cpsi_1eorbit + external linesearch + + real*8 control_time_step,control_fake_mass,ion_disp_energy + external control_time_step,control_fake_mass,ion_disp_energy + integer control_lmbfgs_size + external control_lmbfgs_size + logical control_precondition,band_hfx,ion_disp_on + external control_precondition,band_hfx,ion_disp_on + + real*8 c_geodesic_ABC_energy,c_geodesic_ABC_denergy + external c_geodesic_ABC_energy,c_geodesic_ABC_denergy + + integer cpsi_iptr_psi,c_geodesic_get_U_tag + external cpsi_iptr_psi,c_geodesic_get_U_tag + + dt = control_time_step() +c dte = dt/dsqrt(control_fake_mass()) +c dte = dsqrt(control_fake_mass())/dt +c dte = 1.0d0 + + if (current_iteration.eq.1) then + precondition = control_precondition() + end if + + call Cram_max_npack(npack1) + ispin = cpsi_ispin() + neq(1) = cpsi_neq(1) + neq(2) = cpsi_neq(2) + neall = neq(1)+neq(2) + nbrillq = cpsi_nbrillq() + + +* **** allocate G0, S0 **** + S0_tag = cpsi_data_push_stack(nbrillq,neall,2*npack1) + G0_tag = cpsi_data_push_stack(nbrillq,neall,2*npack1) + G1_tag = cpsi_data_push_stack(nbrillq,neall,2*npack1) + G2_tag = cpsi_data_push_stack(nbrillq,neall,2*npack1) + S0_shift = cpsi_data_get_allptr(S0_tag) + G0_shift = cpsi_data_get_allptr(G0_tag) + +* **** set ptrs **** + Y = cpsi_iptr_psi(1) + U = c_geodesic_get_U_tag() + HY = G1_tag + HU = G2_tag + +* ***** get the initial direction **** + call cpsi_1get_Tgradient(G0_tag,E0) + call cpsi_1get_remainder0(Ermdr) + + if (current_iteration.eq.1) then + +* ***** use the initial gradient for the direction **** + call band_lmbfgs_init(control_lmbfgs_size(),dbl_mb(G0_shift)) + call BGrsm_ff_Copy_tag(G0_tag,S0_tag) + tmin = 10*deltat_min + else + call band_lmbfgs(tmin,dbl_mb(G0_shift),dbl_mb(S0_shift)) + call Pneb_ff_traceall_tag(0,0,G0_tag,S0_tag,sum1) + if (dabs(sum1).gt.1.0d3) call BGrsm_ff_Copy_tag(G0_tag,S0_tag) + end if + + +* ****************************************** +* **** **** +* **** Start of BFGS iteration loop **** +* **** **** +* ****************************************** + it_in = control_it_in() + tole = control_tole() + tolc = control_tolc() + do it=2,it_in + + +* **** initialize the geoedesic line data structure **** + call c_geodesic_start(S0_tag,max_sigma,dE0) + +* **** Copy Hpsi_k to HY then **** +* **** generate Hpsi using psi_k=U and copy Hpsi_k to HU **** +* **** compute A,B,C, , and matrices **** + call c_geodesic_ABC_start(Ermdr,Y,U,HY,HU) + +* ******* line search ********* + if ((tmin.gt.deltat_min).and.(tmin.lt.1.0d4)) then + deltat = tmin + else + deltat = deltat_min + end if +c deltat = 1.0d0 !** debug ** +c Enew = linesearch(0.0d0,E0,dE0,deltat, +c > cpsi_geodesic_energy, +c > cpsi_geodesic_denergy, +c > 0.50d0,tmin,deltae,2) + Enew = linesearch(0.0d0,E0,dE0,deltat, + > c_geodesic_ABC_energy, + > c_geodesic_ABC_denergy, + > 0.50d0,tmin,deltae,2) + call cpsi_geodesic_final(tmin) + call cpsi_2toelectron() + deltac = c_rho_error() + +* **** exit loop early **** + if ((dabs(deltae).lt.tole).and.(deltac.lt.tolc)) then + if (.not.precondition) go to 30 + precondition = .false. + end if + +* **** make psi1 <--- psi2(tmin) **** + call cpsi_2to1() + +* **** get the new gradient - also updates densities**** + call cpsi_1get_Tgradient(G1_tag,E0) + call cpsi_1get_remainder0(Ermdr) + call band_lmbfgs(tmin,dbl_mb(G0_shift),dbl_mb(S0_shift)) + + call Pneb_ff_traceall_tag(0,0,G0_tag,S0_tag,sum1) + if (dabs(sum1).gt.1.0d3) call BGrsm_ff_Copy_tag(G0_tag,S0_tag) + + end do + +* **** initialize the geoedesic line data structure **** + call c_geodesic_start(S0_tag,max_sigma,dE0) + +* **** Copy Hpsi_k to HY then **** +* **** generate Hpsi using psi_k=U and copy Hpsi_k to HU **** +* **** compute A,B,C, , and matrices **** + call c_geodesic_ABC_start(Ermdr,Y,U,HY,HU) + +* ******* line search ********* + if ((tmin.gt.deltat_min).and.(tmin.lt.1.0d4)) then + deltat = tmin + else + deltat = deltat_min + end if +c Enew = linesearch(0.0d0,E0,dE0,deltat, +c > cpsi_geodesic_energy, +c > cpsi_geodesic_denergy, +c > 0.50d0,tmin,deltae,2) + Enew = linesearch(0.0d0,E0,dE0,deltat, + > c_geodesic_ABC_energy, + > c_geodesic_ABC_denergy, + > 0.50d0,tmin,deltae,2) + call cpsi_geodesic_final(tmin) + call cpsi_2toelectron() + deltac = c_rho_error() + +* **** free memory **** + 30 call cpsi_data_pop_stack(G2_tag) + call cpsi_data_pop_stack(G1_tag) + call cpsi_data_pop_stack(G0_tag) + call cpsi_data_pop_stack(S0_tag) + + call cpsi_2to1() + !call cpsi_1toelectron() + call c_rho_2to1() + call c_dng_2to1() +c call cpsi_check() + + eion = ewald_e() + + eorbit = cpsi_1eorbit() + ehartree = c_dng_1ehartree() + exc = c_rho_1exc() + pxc = c_rho_1pxc() + + E(1) = Enew + eion + E(2) = eorbit + E(3) = ehartree + E(4) = exc + E(5) = eion + E(6) = cpsi_1ke() + E(7) = cpsi_1vl() + E(8) = cpsi_1vnl() + E(9) = 2.0d0*ehartree + E(10) = pxc + E(11) = cpsi_1vnlso() + + +* **** HFX terms **** + if (band_HFX()) then + call c_electron_HFX_energies(ehfx,phfx) + E(26) = ehfx + E(27) = phfx + end if + +* **** Dispersion energy **** + if (ion_disp_on()) then + E(33) = ion_disp_energy() + E(1) = E(1) + E(33) + end if + + + return + end + + diff --git a/src/nwpw/band/minimizer/c_cgminimizeABC.F b/src/nwpw/band/minimizer/c_cgminimizeABC.F index d4b58944f8..86a860e2de 100644 --- a/src/nwpw/band/minimizer/c_cgminimizeABC.F +++ b/src/nwpw/band/minimizer/c_cgminimizeABC.F @@ -168,8 +168,6 @@ c > tole,tmin,deltae,1) if ((dabs(deltae).lt.tole).and.(deltac.lt.tolc)) > go to 30 - - * **** transport the previous search directions **** call cpsi_1geodesic_transport( tmin,H0_tag) @@ -178,7 +176,7 @@ c > tole,tmin,deltae,1) * **** get the new gradient - also updates densities**** call cpsi_1get_Tgradient(G1_tag,E0) - + call cpsi_1get_remainder0(Ermdr) sum0 = sum1 diff --git a/src/nwpw/band/minimizer/c_cgsd_energy.F b/src/nwpw/band/minimizer/c_cgsd_energy.F index 70de7cdbcc..eb39a7fd99 100644 --- a/src/nwpw/band/minimizer/c_cgsd_energy.F +++ b/src/nwpw/band/minimizer/c_cgsd_energy.F @@ -194,6 +194,9 @@ c end if else if (minimizer.eq.9) then bfgscount = bfgscount + 1 call c_cgminimizeABC(E,deltae,deltac,bfgscount,it_in) + else if (minimizer.eq.10) then + bfgscount = bfgscount + 1 + call c_bfgsminimizeABC(E,deltae,deltac,bfgscount,minimizer) end if if ((dabs(deltae).gt.dabs(deltae_old)).or. diff --git a/src/nwpw/nwpwlib/control/control.F b/src/nwpw/nwpwlib/control/control.F index 84755284df..8f43a3ffaf 100644 --- a/src/nwpw/nwpwlib/control/control.F +++ b/src/nwpw/nwpwlib/control/control.F @@ -277,6 +277,11 @@ c call nwpwxc_print() > spin_orbit=.false. if (spin_orbit) ispin=2 + if (.not.btdb_get(rtdb,'cpsd:fei',mt_log,1,fei)) + > fei = .false. + if(.not.btdb_cget(rtdb,'cpsd:xyz_filename',1,xyz_filename)) + > call util_file_prefix('xyz',xyz_filename) + * **************************************** * **** cpmd and band_cpmd code: stuff **** diff --git a/src/nwpw/nwpwlib/ion/xyz.F b/src/nwpw/nwpwlib/ion/xyz.F index 185ef314a4..b773ecdcd1 100644 --- a/src/nwpw/nwpwlib/ion/xyz.F +++ b/src/nwpw/nwpwlib/ion/xyz.F @@ -150,5 +150,47 @@ return end +* *************************** +* * * +* * xyz_novelecity_write * +* * * +* *************************** + subroutine xyz_novelocity_write() + implicit none + + real*8 AACONV + parameter (AACONV=0.529177d0) + integer MASTER + parameter (MASTER=0) + + integer taskid,i + +* **** external functions *** + character*2 ion_aname_nocarat + integer ion_nion + real*8 ion_rion, ion_vion + external ion_aname_nocarat + external ion_nion + external ion_rion, ion_vion + integer ierr + + call Parallel_taskid(taskid) + + if (taskid.eq.MASTER) then + write(18,110) ion_nion() + do i=1,ion_nion() + write(18,111) ion_aname_nocarat(i), + > ion_rion(1,i)*AACONV, + > ion_rion(2,i)*AACONV, + > ion_rion(3,i)*AACONV + end do + call util_flush(18) + end if + 110 format(I5/) + 111 format(A2,5x, 6e12.4) + + return + end + diff --git a/src/nwpw/pspw/cpsd/cpsd_input.F b/src/nwpw/pspw/cpsd/cpsd_input.F index 3fd84a7397..7ac330d698 100644 --- a/src/nwpw/pspw/cpsd/cpsd_input.F +++ b/src/nwpw/pspw/cpsd/cpsd_input.F @@ -13,7 +13,7 @@ c c integer ind ! Index of matched directive integer num_dirs ! No. of known directives - parameter (num_dirs = 34) + parameter (num_dirs = 38) character*30 dirs(num_dirs) character*255 test, id @@ -51,11 +51,17 @@ c > 'mult', > 'dft', > 'odft', + > 'fei', + > 'xyz_filename', + > 'ion_motion_filename', + > 'emotion_filename', > 'end'/ character*50 cell_name character*50 input_wavefunction_filename character*50 output_wavefunction_filename + character*50 ion_motion_filename,xyz_filename + character*50 emotion_filename character*50 exchange_correlation logical geometry_optimize,frac_coord,mulliken double precision fake_mass,time_step,rcut @@ -83,7 +89,7 @@ c call cpsd_input_default(rtdb) goto ( 100,100, 200, 300,300, 400,400, 500,500, 600,600, > 700,700, 800,800, 900,900, 1000,1000, 1100,1100, > 1200,1200, 1300,1300, 1400,1400,1400, - > 1500,1600,1700,1800,1900, + > 1500,1600,1700,1800,1900,2000,2100,2200,2300, > 9999) ind call errquit( > 'cpsd_input: unimplemented directive', ind, INPUT_ERR) @@ -329,6 +335,67 @@ c > call errquit('cpsd_input: writing ispin', 0, RTDB_ERR) goto 10 +c +c Fei +c + 2000 if (inp_a(ion_motion_filename)) then + ind = index(ion_motion_filename,' ') - 1 + value = rtdb_cput(rtdb,'cpsd:fei_filename', + > 1,ion_motion_filename(1:ind)) + value = rtdb_cput(rtdb,'nwpw:fei_filename', + > 1,ion_motion_filename(1:ind)) + if (.not.value) + > call errquit('cpsd_input: writing ', 2000, RTDB_ERR) + end if + value = .true. + if (.not.rtdb_put(rtdb,'cpsd:fei',mt_log,1,value)) + > call errquit('cpsd_input: writing fei ', 2000, RTDB_ERR) + if (.not.rtdb_put(rtdb,'nwpw:fei',mt_log,1,value)) + > call errquit('cpsd_input: writing fei ', 2000, RTDB_ERR) + goto 10 + +c +c xyz_filename +c + 2100 if (.not. inp_a(xyz_filename)) + > call errquit( + > 'cpsd_input: failed to read xyz_filename', 0, INPUT_ERR) + ind = index(xyz_filename,' ') - 1 + value = rtdb_cput(rtdb,'cpsd:xyz_filename', + > 1,xyz_filename(1:ind)) + if (.not.value) + > call errquit('cpsd_input: writing', 2100, RTDB_ERR) + goto 10 + +c +c ion_motion_filename +c + 2200 if (.not. inp_a(ion_motion_filename)) + > call errquit( + > 'cpsd_input: failed to read ion_motion_filename', 0, + & INPUT_ERR) + ind = index(ion_motion_filename,' ') - 1 + value = rtdb_cput(rtdb,'cpsd:ion_motion_filename', + > 1,ion_motion_filename(1:ind)) + if (.not.value) + > call errquit('cpsd_input: writing ', 2400, RTDB_ERR) + goto 10 + +c +c emotion_filename +c + 2300 if (.not. inp_a(emotion_filename)) + > call errquit( + > 'cpsd_input: failed to read emotion_filename', 0, + & INPUT_ERR) + ind = index(emotion_filename,' ') - 1 + value = rtdb_cput(rtdb,'cpsd:emotion_filename', + > 1,emotion_filename(1:ind)) + if (.not.value) + > call errquit('cpsd_input: writing ', 2300, RTDB_ERR) + goto 10 + + * ***** add wavefunction to rtdb **** diff --git a/src/nwpw/pspw/cpsd/cpsdv5.F b/src/nwpw/pspw/cpsd/cpsdv5.F index 317710b50d..911c82ac20 100644 --- a/src/nwpw/pspw/cpsd/cpsdv5.F +++ b/src/nwpw/pspw/cpsd/cpsdv5.F @@ -63,7 +63,7 @@ cccc#include "frac_occ.fh" * ***** energy variables **** - real*8 E(60),Egas + real*8 E(60),Egas,dt,emotion_time_shift * real*8 eig(2*nemax) * real*8 hml(2*nemax*nemax) @@ -82,7 +82,7 @@ cccc#include "frac_occ.fh" integer ierr * **** local variables **** - logical oprint,lprint,hprint,mparallelized + logical oprint,lprint,hprint,mparallelized,found,found_bak integer ms,mapping,mapping1d real*8 deltae,deltac,deltar real*8 gx,gy,gz,cx,cy,cz,sum1,sum2 @@ -93,7 +93,7 @@ cccc#include "frac_occ.fh" real*8 w,sumall,virial,dv integer nfft3 parameter (nfft3=32) - character*255 full_filename + character*255 full_filename,full_bak integer tmp1(2) logical value,psi_nogrid,field_exist @@ -106,8 +106,6 @@ cccc#include "frac_occ.fh" logical wvfnc_expander external wvfnc_expander - - * **** external functions **** real*8 psp_zv,psp_rc,ewald_rcut,ion_amass,nwpw_cosmo_qc0 @@ -226,6 +224,7 @@ c#include "Parallel.fh" value = .true. pi = 4.0d0*datan(1.0d0) + dt = 1.0d0 call nwpw_timing_init() call dcopy(60,0.0d0,0,E,1) @@ -745,6 +744,58 @@ c write(luout,1151) (psp_rc(i,ia),i=0,psp_lmax(ia)) ****************** start iterations ********************** * |***************************| +* **** open xyz and MOTION and fei file **** + call xyz_init() ! unit=18 + call MOTION_init(rtdb) ! unit=19 + call fei_init(rtdb) + +* **** open EMOTION file **** + if (.not.btdb_cget(rtdb,'cpsd:emotion_filename',1,filename)) + > call util_file_prefix('emotion',filename) + call util_file_name_noprefix(filename,.false., + > .false., + > full_filename) + if (taskid.eq.MASTER) then + +* **** check for backup file **** + call util_file_name_noprefix('EMOTION99-bak',.false., + > .false., + > full_bak) + inquire(file=full_bak,exist=found_bak) + if (found_bak) then + write(*,*) + write(*,*) "EMOTION99-bak exists:" + i=index(full_bak,' ') + j=index(full_filename,' ') + write(*,*) " Copying ",full_bak(1:i), + > " to ",full_filename(1:j) + write(*,*) + call util_file_copy(full_bak,full_filename) + end if + + emotion_time_shift = 0.0d0 + inquire(file=full_filename,exist=found) + if (found) then + +* **** make a new backup file *** + call util_file_copy(full_filename,full_bak) + + open(unit=31,file=full_filename,form='formatted', + > status='old') + do while (found) + read(31,*,end=100) emotion_time_shift,w,sumall,gx,gy,gz + end do + 100 continue +#if defined(FUJITSU_SOLARIS) || defined(PSCALE) || defined(__crayx1) || defined(GCC46) + backspace 31 +#endif + else + open(unit=31,file=full_filename,form='formatted', + > status='new') + end if + end if + + if (oprint) call current_second(cpu2) if (oprint) call nwpw_message(14) @@ -786,6 +837,16 @@ c write(luout,1151) (psp_rc(i,ia),i=0,psp_lmax(ia)) go to 2 end if + + !**** emotion write **** + write(31,1311) icount*it_in*dt + emotion_time_shift, + > e(1),deltae,deltac,deltar,e(2),e(3),e(4), + > e(5),e(6),e(7),e(8),e(9),e(10) + call util_flush(31) + + !**** xyz write **** + call xyz_novelocity_write() + * ***** define fractional occupation **** if ((icount.lt.it_out).and.fractional.and. > (smeartype.ge.0)) then @@ -822,6 +883,20 @@ c end do if (oprint) CALL nwpw_message(3) if (oprint) call current_second(cpu3) + if (taskid.eq.MASTER) then + close(unit=31) +* **** remove EMOTION backup file *** + call util_file_name_noprefix('EMOTION99-bak',.false., + > .false., + > full_bak) + call util_file_unlink(full_bak) + end if + +* *** close xyz and fei io **** + call xyz_end() + call MOTION_end() + call fei_end() + * |****************************************| *********** produce CHECK file and diagonalize hml ***************** * |****************************************| @@ -1459,6 +1534,7 @@ c call frac_occ_end() 1305 FORMAT(10X, > '==================== iteration =========================') 1310 FORMAT(I8,E20.10,3E15.5) + 1311 format(100e19.10) 1320 FORMAT(' number of electrons: spin up=',F11.5,' down=',F11.5,A) 1321 FORMAT(' total charge of system:',F11.5,A) 1322 FORMAT(' plane-wave part: ',F11.5,' ',F11.5,A) diff --git a/src/nwpw/pspw/cpsd/fei_output.F b/src/nwpw/pspw/cpsd/fei_output.F index 8f3aaa06b4..da804ba9ac 100644 --- a/src/nwpw/pspw/cpsd/fei_output.F +++ b/src/nwpw/pspw/cpsd/fei_output.F @@ -25,13 +25,19 @@ * **** external functions **** logical control_Fei external control_Fei + integer control_code + external control_code if (control_Fei()) then - if (.not.btdb_cget(rtdb,'cpmd:fei_filename',1,filename)) + if ((control_code().eq.1).or.(control_code().eq.13)) then + if (.not.btdb_cget(rtdb,'cpsd:fei_filename',1,filename)) > call util_file_prefix('fei',filename) - + else + if (.not.btdb_cget(rtdb,'cpmd:fei_filename',1,filename)) + > call util_file_prefix('fei',filename) + end if call util_file_name_noprefix(filename,.false., > .false., > full_filename) diff --git a/src/nwpw/pspw/cpsd/inner_loop.F b/src/nwpw/pspw/cpsd/inner_loop.F index 7d14a76532..b2ad3fb8d5 100644 --- a/src/nwpw/pspw/cpsd/inner_loop.F +++ b/src/nwpw/pspw/cpsd/inner_loop.F @@ -94,13 +94,13 @@ c real*8 fion(3,natmx) integer nion1 * ***** external functions **** - logical control_move,psp_semicore,psp_pawexist + logical control_move,psp_semicore,psp_pawexist,control_Fei logical pspw_SIC,pspw_SIC_relaxed logical pspw_HFX,pspw_HFX_relaxed integer ion_nion,control_gga real*8 control_time_step,control_fake_mass,ion_dti real*8 lattice_omega,coulomb_e,ewald_e - external control_move,psp_semicore,psp_pawexist + external control_move,psp_semicore,psp_pawexist,control_Fei external pspw_SIC,pspw_SIC_relaxed external pspw_HFX,pspw_HFX_relaxed external ion_nion,control_gga @@ -206,6 +206,7 @@ c real*8 fion(3,natmx) cosmo1_on = nwpw_cosmo1_on() V_APC_on = pspw_V_APC_on() + call nwpw_timing_start(12) * **** allocate MA local variables **** value = Dneall_m_push_get_block(1,8,tmp_L) @@ -978,6 +979,8 @@ c > - ehartree_atom - ecmp_cmp - ecmp_pw E(1) = E(1) + E(33) end if +* **** output Forces for Fei *** + if (control_Fei()) call fei_output(E(1),dbl_mb(fion(1))) * **** set convergence variables **** diff --git a/src/nwpw/pspw/dplot/ELF.F b/src/nwpw/pspw/dplot/ELF.F index f477b7cd8b..8b43ff8d90 100644 --- a/src/nwpw/pspw/dplot/ELF.F +++ b/src/nwpw/pspw/dplot/ELF.F @@ -66,20 +66,26 @@ > dbl_mb(dns(1)), > dbl_mb(dnsp(1))) do i=1,n2ft3d - dbl_mb(dns(1)+i-1) = log(dn(i)) + dbl_mb(dns(1)+i-1) = log(dn(i)) + write(*,*) "i,dns=",i,dbl_mb(dns(1)+i-1) end do + call D3dB_r_Zero_Ends(1,dbl_mb(dns(1))) * **** elf <-- elf + |del*rho|^2/rho **** - do i=1,n2ft3d - dbl_mb(dns(1)+i-1) = log(dn(i)) - end do + !do i=1,n2ft3d + ! dbl_mb(dns(1)+i-1) = log(dn(i)) + !end do !call D3dB_r_SMul(1,scal1,dbl_mb(dns(1)),dcpl_mb(dng(1))) - call D3dB_r_SMul(1,scal1,dbl_mb(int(dn(1))),dcpl_mb(dng(1))) + !call D3dB_r_SMul(1,scal1,dbl_mb(int(dn(1))),dcpl_mb(dng(1))) + + call D3dB_r_SMul(1,scal1,dn,dcpl_mb(dng(1))) call D3dB_r_SMul(1,scal1,dbl_mb(dns(1)),dcpl_mb(dng2(1))) call D3dB_r_Zero_Ends(1,dcpl_mb(dng(1))) call D3dB_r_Zero_Ends(1,dcpl_mb(dng2(1))) + call D3dB_rc_fft3f(1,dcpl_mb(dng(1))) call D3dB_rc_fft3f(1,dcpl_mb(dng2(1))) + call Pack_c_pack(0,dcpl_mb(dng(1))) call Pack_c_pack(0,dcpl_mb(dng2(1))) call gradient_energy_density2(0,npack0, diff --git a/src/nwpw/pspw/lib/psi/psi.F b/src/nwpw/pspw/lib/psi/psi.F index 27fe05275c..7b1c963498 100644 --- a/src/nwpw/pspw/lib/psi/psi.F +++ b/src/nwpw/pspw/lib/psi/psi.F @@ -2435,6 +2435,7 @@ c + * *********************************** * * * * * psi_3x3_virtual_gradient * @@ -2492,9 +2493,11 @@ c #include "stdio.fh" * **** local variables **** + integer MASTER,taskid + parameter (MASTER=0) logical oprint integer nx,ny,nz,n2ft3d - real*8 Sgg,Sge,See,c1,c2,c3,A,B,C,ttt + real*8 Sgg,Sge,See,c1,c2,c3,A,B,C real*8 H1gg,H1ge,H1gm,H1eg,H1ee,H1em,H1mg,H1me,H1mm real*8 H2gg,H2ge,H2gm,H2eg,H2ee,H2em,H2mg,H2me,H2mm real*8 Hgg,Hge,Hgm,Heg,Hee,Hem,Hmg,Hme,Hmm @@ -2512,9 +2515,14 @@ c real*8 coulomb_e,ewald_e,ion_ion_e external coulomb_e,ewald_e,ion_ion_e +c call Parallel_taskid(taskid) !write(*,*) "ENTER 3x3 GRADIENT" c oprint = .true. +c oprint = (taskid.eq.MASTER) oprint = .false. +c if (oprint) then +c write(*,*) "here here taskid=",taskid +c end if !*** generate psie_r **** call Pack_c_Copy(1,psie,psie_r) @@ -2618,7 +2626,7 @@ c oprint = .true. call coulomb_v(rhogg,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcgg) + call D3dB_r_Copy(1,vcc,vcgg) call D3dB_r_nZero(1,neq(1)+neq(2),hpsig_r) call pspw_potential_HFX(ispin,psig_r,hpsig_r) call pspw_energy_HFX(ispin,psig_r,ehfx,phfx) @@ -2633,19 +2641,13 @@ c oprint = .true. H2ge = -ehfx * *** generate *** - call pspw_et_gen_rho12(ispin,neq,n2ft3d,psig_r,psie_r,rhoge) call D3dB_r_SMul1(1,scal1,rhoge) call D3dB_rc_fft3f(1,rhoge) call Pack_c_pack(0,rhoge) call coulomb_v(rhoge,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcge) - - call pspw_et_gen_rho12(ispin,neq,n2ft3d,psig_r,psie_r,rhoge) - call D3dB_r_SMul1(1,scal1,rhoge) - call D3dB_rc_fft3f(1,rhoge) - call Pack_c_pack(0,rhoge) + call D3dB_r_Copy(1,vcc,vcge) call coulomb_screened_v(rhoge,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) @@ -2662,30 +2664,13 @@ c oprint = .true. call coulomb_v(rhoee,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcee) + call D3dB_r_Copy(1,vcc,vcee) call D3dB_r_nZero(1,neq(1)+neq(2),hpsig_r) call pspw_potential_HFX(ispin,psie_r,hpsig_r) call pspw_energy_HFX(ispin,psie_r,ehfx,phfx) H2ee = H2ee + ehfx * *** generate *** - call pspw_et_gen_rho12(ispin,neq,n2ft3d,psig_r,psie_r,rhoge) - call D3dB_r_SMul1(1,scal1,rhoge) - call D3dB_rc_fft3f(1,rhoge) - call Pack_c_pack(0,rhoge) - call coulomb_v(rhoge,vcc) - call Pack_c_unpack(0,vcc) - call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcge) - - call pspw_et_gen_rho12(ispin,neq,n2ft3d,psig_r,psie_r,rhoge) - call D3dB_r_SMul1(1,scal1,rhoge) - call D3dB_rc_fft3f(1,rhoge) - call Pack_c_pack(0,rhoge) - call coulomb_screened_v(rhoge,vcc) - call Pack_c_unpack(0,vcc) - call D3dB_cr_fft3b(1,vcc) - call D3dB_rr_daxpy(1,-0.5d0,vcc,vcge) call pspw_et_gen_rho(ispin,neq,n2ft3d,psie_r,rhoee) call D3dB_rr_dot(1,rhoee,vcge,H2em) H2em = dsqrt(2.0d0)*(0.50d0*H2em*dv) @@ -2697,7 +2682,7 @@ c oprint = .true. call coulomb_v(rhoee,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcee) + call D3dB_r_Copy(1,vcc,vcee) call coulomb_screened_v(rhoee,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) @@ -2789,15 +2774,15 @@ c Hme = H1me + H2me call coulomb_v(rhoee,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcee) - + call D3dB_r_Copy(1,vcc,vcee) + call D3dB_r_Copy(1,vcee,vc) call D3dB_rr_Mul2(1,psie_r,vc) call D3dB_rr_Sum2(1,hpsig_r,vc) call D3dB_r_SMul1(1,(-4.0d0*c2*c2),vc) - + call D3dB_rr_daxpy(1,(4.0d0*c1*c2),hpsie_r,vc) - + call pspw_et_gen_rho(ispin,neq,n2ft3d,psig_r,rhogg) call D3dB_r_SMul1(1,scal1,rhogg) call D3dB_rc_fft3f(1,rhogg) @@ -2805,7 +2790,7 @@ c Hme = H1me + H2me call coulomb_v(rhogg,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcgg) + call D3dB_r_Copy(1,vcc,vcgg) call coulomb_screened_v(rhogg,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) @@ -2818,7 +2803,7 @@ c Hme = H1me + H2me call coulomb_v(rhoee,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcee) + call D3dB_r_Copy(1,vcc,vcee) call coulomb_screened_v(rhoee,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) @@ -2831,7 +2816,7 @@ c Hme = H1me + H2me call coulomb_v(rhoge,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcge) + call D3dB_r_Copy(1,vcc,vcge) call coulomb_screened_v(rhoge,vcc) call Pack_c_unpack(0,vcc) call D3dB_cr_fft3b(1,vcc) @@ -2840,35 +2825,23 @@ c Hme = H1me + H2me call D3dB_r_Copy(1,vcee,vctmp) call D3dB_rr_Mul2(1,psig_r,vctmp) call D3dB_rr_daxpy(1,(-dsqrt(2.0d0)*2.0d0*c2*c3),vctmp,vc) - + call D3dB_r_Copy(1,vcge,vctmp) call D3dB_rr_Mul2(1,psie_r,vctmp) call D3dB_rr_daxpy(1,(-dsqrt(2.0d0)*2.0d0*c2*c3),vctmp,vc) - +c call D3dB_r_Copy(1,vcgg,vctmp) call D3dB_rr_Mul2(1,psig_r,vctmp) call D3dB_rr_daxpy(1,(-dsqrt(2.0d0)*2.0d0*c3*c1),vctmp,vc) - +c call D3dB_r_Copy(1,vcge,vctmp) call D3dB_rr_Mul2(1,psie_r,vctmp) call D3dB_rr_daxpy(1,(-dsqrt(2.0d0)*2.0d0*c3*c2),vctmp,vc) - call pspw_et_gen_rho(ispin,neq,n2ft3d,psig_r,rhogg) - call D3dB_r_SMul1(1,scal1,rhogg) - call D3dB_rc_fft3f(1,rhogg) - call Pack_c_pack(0,rhogg) - call coulomb_v(rhogg,vcc) - call Pack_c_unpack(0,vcc) - call D3dB_cr_fft3b(1,vcc) - call D3dB_c_Copy(1,vcc,vcgg) - call coulomb_screened_v(rhogg,vcc) - call Pack_c_unpack(0,vcc) - call D3dB_cr_fft3b(1,vcc) - call D3dB_rr_daxpy(1,-0.5d0,vcc,vcgg) call D3dB_r_Copy(1,vcgg,vctmp) call D3dB_rr_Mul2(1,psie_r,vctmp) call D3dB_rr_daxpy(1,(-2.0d0*c3*c3),vctmp,vc) - + call D3dB_rr_daxpy(1,(2.0d0*c3*c3),hpsie_r,vc) end if @@ -2881,13 +2854,6 @@ c Hme = H1me + H2me call Pack_cc_daxpy(1,2.0d0*dsqrt(2.0d0)*c3*c1,h1psig,dEdpsie) call Pack_cc_daxpy(1,2.0d0*c3*c3,h1psie,dEdpsie) -c call Pack_cc_dot(1,psie,dEdpsie,c1) -c write(*,*) "check last gradient=",c1 -c call Pack_cc_dot(1,psig,dEdpsie,c1) -c write(*,*) "check last gradient=",c1 -c call Pack_cc_dot(1,dEdpsie,dEdpsie,c1) -c write(*,*) "check last gradient=",c1 - E = Eci(1) !if (oprint) write(*,*) "c1,c2=",c1,c2 @@ -2895,6 +2861,8 @@ c write(*,*) "check last gradient=",c1 end + + * *********************************** * * * * * psi_3x3ne_virtual_gradient * diff --git a/src/nwpw/pspw/lib/psp/MOTION.F b/src/nwpw/pspw/lib/psp/MOTION.F index 60664e0406..3f7297772f 100644 --- a/src/nwpw/pspw/lib/psp/MOTION.F +++ b/src/nwpw/pspw/lib/psp/MOTION.F @@ -44,18 +44,25 @@ * **** external functions *** real*8 lattice_omega - integer ion_nion,control_it_out + integer ion_nion,control_it_out,control_code external lattice_omega - external ion_nion,control_it_out + external ion_nion,control_it_out,control_code call Parallel_taskid(taskid) - if (.not.btdb_cget(rtdb,'cpmd:ion_motion_filename',1,filename)) - > call util_file_prefix('ion_motion',filename) + if ((control_code().eq.1).or.(control_code().eq.13)) then + if (.not.btdb_cget(rtdb,'cpsd:ion_motion_filename',1,filename)) + > call util_file_prefix('ion_motion',filename) + else + if (.not.btdb_cget(rtdb,'cpmd:ion_motion_filename',1,filename)) + > call util_file_prefix('ion_motion',filename) + end if + call util_file_name_noprefix(filename,.false., > .false., - > full_filename) + > full_filename) + * **** produce MOTION FILE **** if (taskid.eq.MASTER) then diff --git a/src/nwpw/pspw/lib/psp/psp.F b/src/nwpw/pspw/lib/psp/psp.F index b919523c2d..a3f340531b 100644 --- a/src/nwpw/pspw/lib/psp/psp.F +++ b/src/nwpw/pspw/lib/psp/psp.F @@ -266,7 +266,7 @@ c > 'Gijl',Gijl(2),Gijl(1)) do l=0,psp_fmm_lmax do m=-l,l call Tesseral3_rgrid_lm_rl(l,m, - > n2ft3d,dbl_mb(rgrid(1)), + > n2ft3d,dbl_mb(tmp_rgrid(1)), > dbl_mb(psp_fmm_rTlm(1)+lm*n2ft3d)) lm = lm + 1 end do