comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-29 06:35:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Merge pull request #488 from nwchemgit/eric-features

Browse source 
Eric features
This commit is contained in:
Eric Bylaska 2022-02-07 15:40:48 -08:00 • committed by GitHub
commit 5f962f7e4b
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
14 changed files with 545 additions and 83 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/nwpw/band/minimizer/GNUmakefile
View file

@ -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 \

280
src/nwpw/band/minimizer/c_bfgsminimizeABC.F Normal file
View file

@ -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, <S0|H|S0>, <U|H|S0> and <U|H|U> 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, <S0|H|S0>, <U|H|S0> and <U|H|U> 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

4
src/nwpw/band/minimizer/c_cgminimizeABC.F
View file

@ -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

3
src/nwpw/band/minimizer/c_cgsd_energy.F
View file

@ -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.

5
src/nwpw/nwpwlib/control/control.F
View file

@ -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 ****

42
src/nwpw/nwpwlib/ion/xyz.F
View file

@ -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

71
src/nwpw/pspw/cpsd/cpsd_input.F
View file

@ -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 ****

86
src/nwpw/pspw/cpsd/cpsdv5.F
View file

@ -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)

10
src/nwpw/pspw/cpsd/fei_output.F
View file

@ -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)

7
src/nwpw/pspw/cpsd/inner_loop.F
View file

@ -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 ****

16
src/nwpw/pspw/dplot/ELF.F
View file

@ -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,

84
src/nwpw/pspw/lib/psi/psi.F
View file

@ -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 <g|H2|m> ***
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 <e|H2|m> ***
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 *

17
src/nwpw/pspw/lib/psp/MOTION.F
View file

@ -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

2
src/nwpw/pspw/lib/psp/psp.F
View file

@ -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