diff --git a/build_things.sh b/build_things.sh index 83800fa9..f17846c8 100755 --- a/build_things.sh +++ b/build_things.sh @@ -172,6 +172,7 @@ refine_structure thermal_conductivity thermal_conductivity_2023 anharmonic_free_energy +effective_hamiltonian " # only when we have cgal diff --git a/src/anharmonic_free_energy/Makefile b/src/anharmonic_free_energy/Makefile index 2eeeeee8..079460ad 100644 --- a/src/anharmonic_free_energy/Makefile +++ b/src/anharmonic_free_energy/Makefile @@ -3,7 +3,7 @@ CODE=anharmonic_free_energy PROG=../../build/$(CODE)/$(CODE) OBJECT_PATH=../../build/$(CODE)/ -OBJS = $(OBJECT_PATH)main.o $(OBJECT_PATH)options.o $(OBJECT_PATH)energy.o $(OBJECT_PATH)epot.o +OBJS = $(OBJECT_PATH)main.o $(OBJECT_PATH)options.o $(OBJECT_PATH)energy.o LPATH = -L../../lib $(blaslapackLPATH) $(incLPATHhdf) $(incLPATHmpi) IPATH = -I../../inc/libolle -I../../inc/libflap $(blaslapackIPATH) $(incIPATHhdf) $(incIPATHmpi) @@ -24,12 +24,11 @@ clean: .f90.o: $(F90) $(F90FLAGS) -c $< $(LIBS) -$(OBJECT_PATH)main.o: main.f90 $(OBJECT_PATH)options.o $(OBJECT_PATH)energy.o $(OBJECT_PATH)epot.o +$(OBJECT_PATH)main.o: main.f90 $(OBJECT_PATH)options.o $(OBJECT_PATH)energy.o $(F90) $(F90FLAGS) -c main.f90 $(LIBS) -o $@ $(OBJECT_PATH)energy.o: $(OBJECT_PATH)options.o $(F90) $(F90FLAGS) -c energy.f90 $(LIBS) -o $@ $(OBJECT_PATH)options.o: $(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@ -$(OBJECT_PATH)epot.o: epot.f90 - $(F90) $(F90FLAGS) -c epot.f90 $(LIBS) -o $@ + diff --git a/src/anharmonic_free_energy/epot.f90 b/src/anharmonic_free_energy/epot.f90 deleted file mode 100644 index 9085ea65..00000000 --- a/src/anharmonic_free_energy/epot.f90 +++ /dev/null @@ -1,534 +0,0 @@ -module epot -!! Deal with many kinds of potential energy differences -use konstanter, only: r8, lo_pi, lo_twopi, lo_tol, lo_sqtol, lo_status, lo_Hartree_to_eV, lo_kb_hartree -use gottochblandat, only: tochar, walltime, lo_chop, lo_trueNtimes, lo_progressbar_init, & - lo_progressbar, lo_frobnorm, open_file, lo_flattentensor, lo_sqnorm, lo_outerproduct, lo_mean, & - lo_points_on_sphere, lo_mean, lo_stddev -use mpi_wrappers, only: lo_mpi_helper, lo_stop_gracefully -use lo_memtracker, only: lo_mem_helper -!use geometryfunctions, only: lo_inscribed_sphere_in_box -use type_crystalstructure, only: lo_crystalstructure -use type_forceconstant_secondorder, only: lo_forceconstant_secondorder -use type_forceconstant_thirdorder, only: lo_forceconstant_thirdorder -use type_forceconstant_fourthorder, only: lo_forceconstant_fourthorder -use type_mdsim, only: lo_mdsim -implicit none - -private -public :: lo_energy_differences - -type lo_energy_differences - ! forceconstants needed for evaluation - type(lo_forceconstant_secondorder) :: fc2 - type(lo_forceconstant_thirdorder) :: fc3 - type(lo_forceconstant_fourthorder) :: fc4 - real(r8), dimension(:, :, :, :), allocatable :: fcp - logical :: polar = .false. - -contains - procedure :: setup => setup_potential_energy_differences - procedure :: statistical_sampling - procedure :: energies_and_forces -end type - -contains - -!> statistically sample -subroutine statistical_sampling(pot, uc, ss, fc, temperature, mw, mem, verbosity) - !> container for potential energy differences - class(lo_energy_differences), intent(inout) :: pot - !> unitcell - type(lo_crystalstructure), intent(in) :: uc - !> supercell - type(lo_crystalstructure), intent(in) :: ss - !> second order forceconstant - type(lo_forceconstant_secondorder), intent(in) :: fc - !> temperature - real(r8), intent(in) :: temperature - !> MPI helper - type(lo_mpi_helper), intent(inout) :: mw - !> memory tracker - type(lo_mem_helper), intent(inout) :: mem - !> talk a lot? - integer, intent(in) :: verbosity - - type(lo_crystalstructure) :: p - integer, parameter :: ninner = 20, nouter = 40 - integer :: initer, outiter, nstep, ctr, i - real(r8), dimension(:, :), allocatable :: ebuf, rbuf - real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp - real(r8) :: e2, e3, e4, ep, deltaE, expdeltaE, ikbt, tomev, emean - - ! Copy of structure to work with - p = ss - - ! Total number of steps we are doing? - nstep = ninner*nouter - - ! Some helper space - call mem%allocate(ebuf, [nstep, 3], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%allocate(rbuf, [3, mw%n], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - ebuf = 0.0_r8 - rbuf = 0.0_r8 - - ! Dummy space for force - call mem%allocate(f2, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%allocate(f3, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%allocate(f4, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%allocate(fp, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - f2 = 0.0_r8 - f3 = 0.0_r8 - f4 = 0.0_r8 - fp = 0.0_r8 - - ! Thermal factor? - ikbt = 1.0_r8/(temperature*lo_kb_hartree) - tomev = lo_Hartree_to_eV*1000 - - ctr = 0 - outloop: do outiter = 1, nouter - inloop: do initer = 1, ninner - ctr = ctr + 1 - - ! Reset structure - p%u = 0.0_r8 - p%v = 0.0_r8 - p%r = ss%r - ! Get new structure - call pot%fc2%initialize_cell(p, uc, fc, temperature, .true., .false., -1.0_r8, mw, nosync=.true.) - ! Calculate the energy - call pot%energies_and_forces(p%u, e2, e3, e4, ep, f2, f3, f4, fp) - e2 = e2 !/ss%na - e3 = e3 !/ss%na - e4 = e4 !/ss%na - ep = ep !/ss%na - - ! H1-H0 difference - deltaE = e3 + e4 - expdeltaE = exp(-deltaE*ikbt) - ebuf(ctr, 1) = deltaE - ebuf(ctr, 2) = expdeltaE - end do inloop - - ! Get new mean energy - emean = lo_mean(ebuf(1:ctr, 1))/mw%n - call mw%allreduce('sum', emean) - do i = 1, ctr - ebuf(i, 3) = ebuf(i, 1) - emean - ebuf(i, 3) = exp(-ebuf(i, 3)*ikbt) - end do - - ! Print a little summary - rbuf = 0.0_r8 - rbuf(1, mw%r + 1) = lo_mean(ebuf(1:ctr, 1))/ss%na - rbuf(2, mw%r + 1) = -log(lo_mean(ebuf(1:ctr, 2)))/ikbt/ss%na - rbuf(3, mw%r + 1) = -log(lo_mean(ebuf(1:ctr, 3)))/ikbt/ss%na + emean/ss%na - call mw%allreduce('sum', rbuf) - - if (verbosity .gt. 0) then - write (*, *) ctr, lo_mean(rbuf(3, :))*tomev, lo_mean(rbuf(2, :))*tomev, lo_stddev(rbuf(1, :))*tomev, lo_stddev(rbuf(2, :))*tomev - end if - end do outloop - - call mem%deallocate(ebuf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%deallocate(rbuf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%deallocate(f2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%deallocate(f3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%deallocate(f4, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - call mem%deallocate(fp, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) -end subroutine - -! !> calculate energies and forces for a given configuration -subroutine energies_and_forces(pot, u, e2, e3, e4, ep, f2, f3, f4, fp) - !> container for potential energy differences - class(lo_energy_differences), intent(in) :: pot - !> displacements - real(r8), dimension(:, :), intent(in) :: u - !> energies - real(r8), intent(out) :: e2, e3, e4, ep - !> forces - real(r8), dimension(:, :), intent(out) :: f2, f3, f4, fp - - real(r8), dimension(3, 3, 3, 3) :: m4 - real(r8), dimension(3, 3, 3) :: m3 - real(r8), dimension(3, 3) :: m2 - real(r8), dimension(3) :: v0, u2, u3, u4 - integer :: a1, a2, a3, a4, i1, i2, i3, i4, i - - e2 = 0.0_r8 - e3 = 0.0_r8 - e4 = 0.0_r8 - ep = 0.0_r8 - f2 = 0.0_r8 - f3 = 0.0_r8 - f4 = 0.0_r8 - fp = 0.0_r8 - - ! pair term - do a1 = 1, size(u, 2) - v0 = 0.0_r8 - do i1 = 1, pot%fc2%atom(a1)%n - a2 = pot%fc2%atom(a1)%pair(i1)%i2 - m2 = pot%fc2%atom(a1)%pair(i1)%m - v0 = v0 - matmul(m2, u(:, a2)) - end do - e2 = e2 - dot_product(u(:, a1), v0)*0.5_r8 - f2(:, a1) = v0 - end do - - ! polar term - if (pot%polar) then - ep = 0.0_r8 - do a1 = 1, size(u, 2) - v0 = 0.0_r8 - do a2 = 1, size(u, 2) - v0 = v0 - matmul(pot%fcp(:, :, a1, a2), u(:, a2)) - end do - ep = ep - dot_product(u(:, a1), v0)*0.5_r8 - fp(:, a1) = v0 - end do - end if - - ! triplet term - if (pot%fc3%na .eq. size(u, 2)) then - do a1 = 1, size(u, 2) - v0 = 0.0_r8 - do i = 1, pot%fc3%atom(a1)%n - m3 = pot%fc3%atom(a1)%triplet(i)%m - a2 = pot%fc3%atom(a1)%triplet(i)%i2 - a3 = pot%fc3%atom(a1)%triplet(i)%i3 - u2 = u(:, a2) - u3 = u(:, a3) - do i1 = 1, 3 - do i2 = 1, 3 - do i3 = 1, 3 - v0(i1) = v0(i1) - m3(i1, i2, i3)*u2(i2)*u3(i3) - end do - end do - end do - end do - v0 = v0*0.5_r8 - f3(:, a1) = v0 - e3 = e3 - dot_product(v0, u(:, a1))/3.0_r8 - end do - end if - - ! quartet term - if (pot%fc4%na .eq. size(u, 2)) then - do a1 = 1, size(u, 2) - v0 = 0.0_r8 - do i = 1, pot%fc4%atom(a1)%n - m4 = pot%fc4%atom(a1)%quartet(i)%m - a2 = pot%fc4%atom(a1)%quartet(i)%i2 - a3 = pot%fc4%atom(a1)%quartet(i)%i3 - a4 = pot%fc4%atom(a1)%quartet(i)%i4 - u2 = u(:, a2) - u3 = u(:, a3) - u4 = u(:, a4) - do i1 = 1, 3 - do i2 = 1, 3 - do i3 = 1, 3 - do i4 = 1, 3 - v0(i1) = v0(i1) - m4(i1, i2, i3, i4)*u2(i2)*u3(i3)*u4(i4) - end do - end do - end do - end do - end do - v0 = v0/6.0_r8 - f4(:, a1) = v0 - e4 = e4 - dot_product(v0, u(:, a1))/4.0_r8 - end do - end if -end subroutine - -!> Calculate potential energy differences in several ways -subroutine setup_potential_energy_differences(pot, uc, ss, fc2, fc3, fc4, mw, verbosity) - !> container for potential energy differences - class(lo_energy_differences), intent(out) :: pot - !> unitcell - type(lo_crystalstructure), intent(inout) :: uc - !> supercell - type(lo_crystalstructure), intent(inout) :: ss - !> second order forceconstant - type(lo_forceconstant_secondorder), intent(in) :: fc2 - !> third order forceconstant - type(lo_forceconstant_thirdorder), intent(in) :: fc3 - !> fourth order forceconstant - type(lo_forceconstant_fourthorder), intent(in) :: fc4 - !> mpi things - type(lo_mpi_helper), intent(inout) :: mw - !> how much to talk - integer, intent(in) :: verbosity - - real(r8) :: timer, t0, t1 - - timer = walltime() - t0 = timer - t1 = timer - - if (verbosity .gt. 0) then - write (*, *) '' - write (*, *) 'PREPARING POTENTIAL ENERGY DIFFERENCES' - end if - - call fc2%remap(uc, ss, pot%fc2) - if (verbosity .gt. 0) then - t1 = walltime() - write (*, *) '... remapped second order (', tochar(t1 - t0), ')' - t0 = t1 - end if - - if (fc3%na .gt. 0) then - call fc3%remap(uc, ss, pot%fc3) - if (verbosity .gt. 0) then - t1 = walltime() - write (*, *) '... remapped third order (', tochar(t1 - t0), ')' - t0 = t1 - end if - end if - - if (fc4%na .gt. 0) then - call fc4%remap(uc, ss, pot%fc4) - if (verbosity .gt. 0) then - t1 = walltime() - write (*, *) '... remapped fourth order (', tochar(t1 - t0), ')' - t0 = t1 - end if - end if - - if (fc2%polar) then - allocate (pot%fcp(3, 3, ss%na, ss%na)) - pot%fcp = 0.0_r8 - call fc2%supercell_longrange_dynamical_matrix_at_gamma(ss, pot%fcp, 1E-15_r8) - pot%polar = .true. - else - pot%polar = .false. - end if - if (verbosity .gt. 0) then - t1 = walltime() - write (*, *) '... built polar forceconstant (', tochar(t1 - t0), ')' - t0 = t1 - end if -end subroutine - -! !> remove longrange interactions -! subroutine potential_energy_difference(sim,uc,ss,fc,fct,fcf,Uref,mw,verbosity) -! !> force-displacement data -! type(lo_mdsim), intent(inout) :: sim -! !> unitcell -! type(lo_crystalstructure), intent(inout) :: uc -! !> supercell -! type(lo_crystalstructure), intent(inout) :: ss -! !> second order forceconstant -! type(lo_forceconstant_secondorder), intent(in) :: fc -! !> third order forceconstant -! type(lo_forceconstant_thirdorder), intent(in) :: fct -! !> fourth order forceconstant -! type(lo_forceconstant_fourthorder), intent(in) :: fcf -! !> reference energy -! real(flyt), intent(in) :: Uref -! !> mpi things -! type(lo_mpi_helper), intent(inout) :: mw -! !> how much to talk -! integer, intent(in) :: verbosity -! -! type(lo_forceconstant_secondorder) :: fcss -! type(lo_forceconstant_thirdorder) :: fctss -! type(lo_forceconstant_fourthorder) :: fcfss -! real(flyt), dimension(:,:,:,:), allocatable :: forceconstant -! real(flyt), dimension(:,:), allocatable :: fp,f2,f3,f4,de -! ! real(flyt), dimension(:,:), allocatable :: f,de -! ! real(flyt), dimension(3,3) :: m0,m1 -! ! real(flyt), dimension(3) :: v0 -! ! real(flyt) :: epot,epot_2,epot_3,epot_4,epot_dd -! ! integer :: a1,a2,t,u -! -! init: block -! real(flyt), dimension(3,3) :: m0,m1 -! real(flyt) :: t0,f0 -! integer :: a1,a2 -! if ( mw%talk ) then -! write(*,*) '' -! write(*,*) 'CALCULATING POTENTIAL ENERGY DIFFERENCES' -! endif -! ! Remap the forceconstants to the supercell -! call ss%classify('supercell',uc) -! call fc%remap(uc,ss,fcss) -! call fct%remap(uc,ss,fctss) -! call fcf%remap(uc,ss,fcfss) -! if ( mw%talk ) write(*,*) '... remapped forceconstants' -! ! Get the dipole-dipole forceconstant -! if ( fc%polar ) then -! lo_allocate(forceconstant(3,3,ss%na,ss%na)) -! call fc%supercell_longrange_dynamical_matrix_at_gamma(ss,forceconstant,1E-15_flyt) -! ! small sanity-check, just in case -! f0=0.0_flyt -! do a1=1,ss%na -! do a2=a1,ss%na -! m0=forceconstant(:,:,a1,a2) -! m1=transpose(forceconstant(:,:,a2,a1)) -! f0=f0+lo_frobnorm(m0-m1) -! enddo -! enddo -! f0=f0/(ss%na**2) -! if ( f0 .gt. lo_tol ) then -! write(*,*) 'ERROR, non-Hermitian electrostatic forceconstant: ',f0 -! write(*,*) 'If this keeps happening I should not hard-code the Ewald tolerance.' -! stop -! endif -! endif -! -! ! some space for energies and forces -! lo_allocate(de(sim%nt,5)) -! lo_allocate(fp(3,ss%na)) -! lo_allocate(f2(3,ss%na)) -! lo_allocate(f3(3,ss%na)) -! lo_allocate(f4(3,ss%na)) -! de=0.0_flyt -! fp=0.0_flyt -! f2=0.0_flyt -! f3=0.0_flyt -! f4=0.0_flyt -! end block init -! -! forceenergy: block -! real(flyt), dimension(3,3,3,3) :: m4 -! real(flyt), dimension(3,3,3) :: m3 -! real(flyt), dimension(3,3) :: m2 -! real(flyt), dimension(3) :: u1,u2,u3,u4,v0 -! real(flyt) :: energy -! integer :: a1,a2,a3,a4,i1,i2,i3,i4,i,j,k,l,t -! -! do t=1,sim%nt -! ! get the raw energy -! de(t,1)=sim%stat%potential_energy(t) -! -! ! secondorder energy, first the polar term (if applicable) -! energy=0.0_flyt -! if ( fc%polar ) then -! do a1=1,ss%na -! v0=0.0_flyt -! do a2=1,ss%na -! v0=v0+matmul(forceconstant(:,:,a2,a1),u(:,a2,t)) -! enddo -! energy=energy+dot_product(u(:,a1,t),v0)*0.5_flyt -! enddo -! endif -! ! then the pair term, always there -! do a1=1,ss%na -! v0=0.0_flyt -! do i1=1,fcss%atom(a1)%n -! a2=fcss%atom(a1)%pair(i1)%i2 -! m2=fcss%atom(a1)%pair(i1)%m -! v0=v0+matmul(m2,u(:,a2,t)) -! enddo -! energy=energy+dot_product(u(:,a1,t),v0)*0.5_flyt -! enddo -! de(t,2)=energy -! -! ! triplet term -! f3=0.0_flyt -! energy=0.0_flyt -! do a1=1,fcss%na -! v0=0.0_flyt -! do i=1,fctss%atom(a1)%n -! m3=fctss%atom(a1)%triplet(i)%m -! a2=fctss%atom(a1)%triplet(i)%i2 -! a3=fctss%atom(a1)%triplet(i)%i3 -! u2=u(:,a2,t) -! u3=u(:,a3,t) -! do i1=1,3 -! do i2=1,3 -! do i3=1,3 -! v0(i1)=v0(i1)-m3(i1,i2,i3)*u2(i2)*u3(i3) -! enddo -! enddo -! enddo -! enddo -! f3(:,a1)=v0*0.5_flyt -! energy=energy-dot_product(f3(:,a1),u(:,a1,t))/3.0_flyt -! enddo -! de(t,3)=energy -! -! ! quartet term -! f4=0.0_flyt -! energy=0.0_flyt -! do a1=1,fcfss%na -! v0=0.0_flyt -! do i=1,fcfss%atom(a1)%n -! m4=fcfss%atom(a1)%quartet(i)%m -! a2=fcfss%atom(a1)%quartet(i)%i2 -! a3=fcfss%atom(a1)%quartet(i)%i3 -! a4=fcfss%atom(a1)%quartet(i)%i4 -! u2=u(:,a2,t) -! u3=u(:,a3,t) -! u4=u(:,a4,t) -! do i1=1,3 -! do i2=1,3 -! do i3=1,3 -! do i4=1,3 -! v0(i1)=v0(i1)-m4(i1,i2,i3,i4)*u2(i2)*u3(i3)*u4(i4) -! enddo -! enddo -! enddo -! enddo -! enddo -! f4(:,a1)=v0/6.0_flyt -! energy=energy-dot_product(f4(:,a1),u(:,a1,t))/4.0_flyt -! enddo -! de(t,4)=energy -! enddo -! end block forceenergy -! -! if ( mw%talk ) then -! write(*,*) ' raw avg potential energy:',lo_mean(de(:,1)) -! write(*,*) ' E-E2:',lo_mean(de(:,1)-de(:,2)) -! write(*,*) ' E-E2-E3:',lo_mean(de(:,1)-de(:,2)-de(:,3)) -! write(*,*) ' E-E2-E3-E4:',lo_mean(de(:,1)-de(:,2)-de(:,3)-de(:,4)) -! endif -! -! -! !! Energy per atom -! !de=1000*de/sim%na -! ! -! !if ( mw%talk ) then -! !write(*,*) 'RAW:',lo_mean(de(:,1)),lo_stddev(de(:,1)) -! !write(*,*) 'RAW:',lo_mean(de(:,1)-de(:,3)),lo_stddev(de(:,1)-de(:,3)) -! !write(*,*) 'RAW:',lo_mean(de(:,1)-de(:,2)-de(:,3)),lo_stddev(de(:,1)-de(:,2)-de(:,3)) -! !write(*,*) 'RAW:',lo_mean(de(:,1)-de(:,2)-de(:,3)-de(:,4)),lo_stddev(de(:,1)-de(:,2)-de(:,3)-de(:,4)) -! !write(*,*) 'RAW:',lo_mean(de(:,1)-de(:,2)-de(:,3)-de(:,4)-de(:,5)),lo_stddev(de(:,1)-de(:,2)-de(:,3)-de(:,4)-de(:,5)) -! !u=open_file('out','outfile.U0') -! ! do t=1,sim%nt -! ! write(u,"(6(1X,E18.12))") de(t,:)/1000*sim%na,(de(t,1)-de(t,2)-de(t,3)-de(t,4)-de(t,5))/1000*sim%na -! ! enddo -! !close(u) -! !endif -! -! ! t0=walltime() -! ! call lo_progressbar_init() -! ! ! Subtract forces -! ! lo_allocate(f(3,ss%na)) -! ! do t=1,sim%nt -! ! ! calculate the long-range forces -! ! f=0.0_flyt -! ! do a1=1,ss%na -! ! do a2=1,ss%na -! ! v0=matmul(forceconstant(:,:,a1,a2),u(:,a2,t)) -! ! f(:,a1)=f(:,a1)-v0 -! ! enddo -! ! enddo -! ! ! sanity check that they add up to zero -! ! if ( abs(sum(f)) .gt. lo_sqtol ) then -! ! write(*,*) '' -! ! write(*,*) 'ERROR: dipole-dipole forces do not add up to zero!' -! ! stop -! ! endif -! ! ! subtract them -! ! sim%f(:,:,t)=sim%f(:,:,t)-f -! ! if( lo_trueNtimes(t,50,sim%nt) ) call lo_progressbar(' ... subtracting dipole-dipole forces',t,sim%nt) -! ! enddo -! ! call lo_progressbar(' ... subtracting dipole-dipole forces',sim%nt,sim%nt,walltime()-t0) -! end subroutine - -end module diff --git a/src/anharmonic_free_energy/main.f90 b/src/anharmonic_free_energy/main.f90 index f815e238..3f57a156 100644 --- a/src/anharmonic_free_energy/main.f90 +++ b/src/anharmonic_free_energy/main.f90 @@ -14,11 +14,11 @@ program anharmonic_free_energy use type_qpointmesh, only: lo_qpoint_mesh, lo_generate_qmesh, lo_read_qmesh_from_file, lo_fft_mesh use type_phonon_dispersions, only: lo_phonon_dispersions use lo_phonon_bandstructure_on_path, only: lo_phonon_bandstructure +use lo_epot, only: lo_energy_differences use type_phonon_dos, only: lo_phonon_dos use options, only: lo_opts use energy, only: perturbative_anharmonic_free_energy -use epot, only: lo_energy_differences implicit none diff --git a/src/effective_hamiltonian/Makefile b/src/effective_hamiltonian/Makefile new file mode 100644 index 00000000..9eb9b4e8 --- /dev/null +++ b/src/effective_hamiltonian/Makefile @@ -0,0 +1,31 @@ +include Makefile.inc +CODE=effective_hamiltonian +PROG=../../build/$(CODE)/$(CODE) +OBJECT_PATH=../../build/$(CODE)/ + +OBJS = $(OBJECT_PATH)main.o $(OBJECT_PATH)options.o + +LPATH = -L../../lib $(blaslapackLPATH) $(incLPATHhdf) $(incLPATHmpi) +IPATH = -I../../inc/libolle -I../../inc/libflap $(blaslapackIPATH) $(incIPATHhdf) $(incIPATHmpi) +LIBS = ../../lib/libolle.a -lflap $(blaslapackLIBS) $(incLIBShdf) $(incLIBSmpi) + +#OPT = -O0 -fbacktrace -fcheck=all -Wall +F90 = $(FC) $(LPATH) $(IPATH) $(MODULE_FLAG) $(OBJECT_PATH) +F90FLAGS = $(OPT) $(MODS) $(LIBS) + +all: $(PROG) + +$(PROG): $(OBJS) + $(F90) $(LDFLAGS) -o $@ $(OBJS) $(LIBS) + +clean: + rm -f $(PROG) $(OBJS) *.mod + +.f90.o: + $(F90) $(F90FLAGS) -c $< $(LIBS) + +$(OBJECT_PATH)main.o: main.f90 $(OBJECT_PATH)options.o + $(F90) $(F90FLAGS) -c main.f90 -o $@ +$(OBJECT_PATH)options.o: + $(F90) $(F90FLAGS) -c options.f90 -o $@ + diff --git a/src/effective_hamiltonian/main.f90 b/src/effective_hamiltonian/main.f90 new file mode 100644 index 00000000..7ff1f8a4 --- /dev/null +++ b/src/effective_hamiltonian/main.f90 @@ -0,0 +1,295 @@ +#include "precompilerdefinitions" +program effective_hamiltonian +!!{!src/effective_hamiltonian/manual.md!} +use konstanter, only: r8, lo_tol, lo_kb_hartree, lo_bohr_to_A, lo_frequency_Hartree_to_THz, lo_eV_to_Hartree, & + lo_exitcode_physical, lo_exitcode_param, lo_temperaturetol, lo_Hartree_to_eV +use mpi_wrappers, only: lo_mpi_helper +use lo_memtracker, only: lo_mem_helper +use gottochblandat, only: tochar, walltime, lo_stop_gracefully, open_file, lo_progressbar_init, & + lo_progressbar, lo_does_file_exist, lo_trueNtimes, lo_mean + +use type_forceconstant_secondorder, only: lo_forceconstant_secondorder +use type_forceconstant_thirdorder, only: lo_forceconstant_thirdorder +use type_forceconstant_fourthorder, only: lo_forceconstant_fourthorder +use type_mdsim, only: lo_mdsim + +use lo_epot, only: lo_energy_differences + +use type_crystalstructure, only: lo_crystalstructure +use options, only: lo_opts + +implicit none +type(lo_opts) :: opts +type(lo_crystalstructure) :: ss, uc +type(lo_energy_differences) :: pot +type(lo_mdsim) :: sim + +type(lo_forceconstant_secondorder) :: fc2 +type(lo_forceconstant_thirdorder) :: fc3 +type(lo_forceconstant_fourthorder) :: fc4 +! POLAR IFCS? + +type(lo_mpi_helper) :: mw +type(lo_mem_helper) :: mem +real(r8), dimension(:, :), allocatable :: pebuf +! real(r8), dimension(:), allocatable :: kebuf + +logical :: generate_configs = .false. + + +call opts%parse() +call mw%init() +call mem%init() + +init: block + + integer :: f, i, j, l, readrank + logical :: readonthisrank, mpiparallel + real(r8) :: t0 + + + if (.not. mw%talk) opts%verbosity = -100 + + if (opts%nconf .gt. 0) generate_configs = .true. + + if (mw%talk) then + write (*, *) 'Recap of the parameters governing the calculation' + write (*, '(1X,A40,L3)') 'Thirdorder contribution ', opts%thirdorder + write (*, '(1X,A40,L3)') 'Fourthorder contribution ', opts%fourthorder + write (*, '(1X,A40,I5)') 'Stride ', opts%stride + write(*, '(1X,A40,L3)') 'Generate canonical configs ', generate_configs + if(generate_configs) then + write(*, '(1X,A40,I8)') 'Num Configs ', opts%nconf + write (*, '(1X,A40,L3)') 'Quantum configurations ', opts%quantum + write (*, '(1X,A40,F20.12)') 'Temperature ', opts%temperature + end if + end if + + ! Read structures + if (mw%talk) write (*, *) '... reading infiles' + call ss%readfromfile('infile.ssposcar', verbosity=opts%verbosity) + call uc%readfromfile('infile.ucposcar', verbosity=opts%verbosity) + + ! Match the supercell to the unitcell, always a good idea + call uc%classify('spacegroup', timereversal=.true.) + call ss%classify('supercell', uc) + + call fc2%readfromfile(uc, 'infile.forceconstant', mem, verbosity=-1) + if (opts%thirdorder) call fc3%readfromfile(uc, 'infile.forceconstant_thirdorder') + if (opts%fourthorder) call fc4%readfromfile(uc, 'infile.forceconstant_fourthorder') + if (mw%talk) write (*, *) '... read forceconstants' + + if (generate_configs) then + if ((opts%quantum .eqv. .false.) .and. (opts%temperature .lt. lo_temperaturetol)) then + call lo_stop_gracefully(['For classical statistics temperature has to be nonzero'], lo_exitcode_physical, __FILE__, __LINE__) + end if + if (mw%talk .and. opts%quantum) write (*, *) '... will generate quantum configurations' + if (mw%talk .and. opts%quantum) write (*, *) '... will generate classical configurations' + + if (opts%temperature .lt. lo_temperaturetol) then + call lo_stop_gracefully(['If nconf is passed, the temperature flag must be set'], lo_exitcode_param, __FILE__, __LINE__) + end if + + end if + + call pot%setup(uc, ss, fc2, fc3, fc4, mw, opts%verbosity + 1) + if (mw%talk) write (*, *) '... setup potential energy calculator' + + if (.not. generate_configs) then + ! If packed simulation is there might as well read it + if (lo_does_file_exist('infile.sim.hdf5')) then + call sim%read_from_hdf5('infile.sim.hdf5', verbosity=opts%verbosity + 2, stride=opts%stride) + else + call sim%read_from_file(verbosity = opts%verbosity + 2, stride = opts%stride, & + magnetic=.false., dielectric=.false., nrand=-1, mw=mw) + end if + if (mw%talk) write (*, *) '... parsed simulation data' + end if + +end block init + + + + +energy : block + + integer :: i, u + real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp + real(r8) :: e2, e3, e4, ep, total_energy, to_ev_per_atom + + + if (generate_configs) then + + if (mw%talk) write (*, *) '... generating canonical configurations' + call mem%allocate(pebuf, [opts%nconf, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + pebuf = 0.0_r8 + + call pot%statistical_sampling(uc, ss, fc2, opts%nconf, opts%temperature, opts%quantum, pebuf, mw, mem, opts%verbosity) + + else + + call mem%allocate(pebuf, [sim%nt, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + pebuf = 0.0_r8 + + ! call mem%allocate(kebuf, [sim%nt], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + ! kebuf = 0.0_r8 + + ! Dummy space for force + call mem%allocate(f2, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(f3, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(f4, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(fp, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + f2 = 0.0_r8 + f3 = 0.0_r8 + f4 = 0.0_r8 + fp = 0.0_r8 + + ! how do I add a progress bar here without race condition? + do i = 1, sim%nt + + if (mod(i, mw%n) .ne. mw%r) cycle + + ! Calculate the energy, e2/e3/e4/ep are zeroed inside of this call + call pot%energies_and_forces(sim%u(:,:,i), e2, e3, e4, ep, f2, f3, f4, fp) + + pebuf(i, 1) = e2 + pebuf(i, 2) = e3 + pebuf(i, 3) = e4 + pebuf(i, 4) = ep + end do + + call mw%allreduce('sum', pebuf) + + end if + + if (mw%talk) write (*, *) '... calculated energies' + +end block energy + +! Modified from anharmonic_free_energy +epotthings: block + real(r8), dimension(3, 5) :: cumulant + real(r8), dimension(:, :), allocatable :: ediff + real(r8) :: inverse_kbt, T_actual, U0, total_energy, hartree_to_mev, to_mev_per_atom + integer :: i, u + + to_mev_per_atom = 1000*lo_Hartree_to_eV / real(ss%na, r8) + + if (generate_configs) then + T_actual = opts%temperature + else + T_actual = sim%temperature_thermostat + end if + + + if (.not. generate_configs) then + ! Calculate the baseline energy + allocate (ediff(sim%nt, 5)) + ediff = 0.0_r8 + + do i = 1, sim%nt + if (mod(i, mw%n) .ne. mw%r) cycle + ediff(i, 1) = sim%stat%potential_energy(i) + ediff(i, 2) = sim%stat%potential_energy(i) - pebuf(i, 1) + ediff(i, 3) = sim%stat%potential_energy(i) - pebuf(i, 1) - pebuf(i, 4) + ediff(i, 4) = sim%stat%potential_energy(i) - pebuf(i, 1) - pebuf(i, 4) - pebuf(i, 2) + ediff(i, 5) = sim%stat%potential_energy(i) - pebuf(i, 1) - pebuf(i, 4) - pebuf(i, 2) - pebuf(i, 3) + end do + call mw%allreduce('sum', ediff) + + if (T_actual .gt. 1E-5_r8) then + inverse_kbt = 1.0_r8/lo_kb_Hartree/T_actual + else + inverse_kbt = 0.0_r8 + end if + + ! Compute the first and second order cumulants + do i = 1, 5 + cumulant(1, i) = lo_mean(ediff(:, i)) + cumulant(2, i) = lo_mean((ediff(:, i) - cumulant(1, i))**2) + cumulant(2, i) = cumulant(2, i)*inverse_kbt*0.5_r8 + cumulant(3, i) = lo_mean((ediff(:, i) - cumulant(1, i))**3) + cumulant(3, i) = cumulant(3, i)*inverse_kbt**2/6.0_r8 + end do + + ! And normalize it to be per atom + cumulant = cumulant*to_mev_per_atom + if (mw%talk) then + u = open_file('out', 'outfile.cumulants') + + write (u, *) '# Temperature (K) : ', T_actual + write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) + write (u, *) '# Potential energy [eV / atom]:' + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E): ', cumulant(1, 1), ' upper bound:', cumulant(2, 1) + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2): ', cumulant(1, 2), ' upper bound:', cumulant(2, 2) + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar): ', cumulant(1, 3), ' upper bound:', cumulant(2, 3) + if (opts%thirdorder) then + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar-E3): ', cumulant(1, 4), ' upper bound:', cumulant(2, 4) + end if + if (opts%fourthorder) then + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar-E3-E4): ', cumulant(1, 5), ' upper bound:', cumulant(2, 5) + end if + + close(u) + write (*, '(A)') ' ... cumulants writen to `outfile.cumulants`' + end if + + if(opts%thirdorder .and. (.not. opts%fourthorder)) then + U0 = cumulant(1,4) + else if(opts%fourthorder) then + U0 = cumulant(1,5) + else + U0 = cumulant(1,3) + end if + + ! Now that we have U0 we can get E_total_tdep + if (mw%talk) then + u = open_file('out', 'outfile.energies') + write (u, '(A,A)') '# Unit: ', 'meV/atom' + write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) + write (u, *) '# U0 [meV/atom]: ', U0 + write (u, "(A)") '# conf Etotal_actual Etotal_tdep Epolar & + &Epair Etriplet Equartet' + + do i = 1, sim%nt + total_energy = sum(pebuf(i,:))*to_mev_per_atom + U0 + write (u, "(1X,I8,6(2X,E20.12))") i, sim%stat%potential_energy(i)*to_mev_per_atom, total_energy, pebuf(i, 4)*to_mev_per_atom, & + pebuf(i, 1)*to_mev_per_atom, pebuf(i, 2)*to_mev_per_atom, & + pebuf(i, 3)*to_mev_per_atom + end do + + ! close outfile.energies + close (u) + + write (*, '(A)') ' ... energies writen to `outfile.energies`' + end if + + ! Don't have sim%stat%potential energy when generating canonical_configs + ! So we cannot get U0 or cumulants + else + if (mw%talk) then + u = open_file('out', 'outfile.energies') + write (u, '(A,A)') '# Unit: ', 'meV/atom' + write (u, *) '# Temperature (K) : ', T_actual + write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) + write (u, "(A)") '# True potential energy unknown when using canonical configs, cannot calculate U0 or other cumulants.' + write (u, "(A)") '# conf Epolar & + &Epair Etriplet Equartet' + + do i = 1, opts%nconf + write (u, "(1X,I8,4(2X,E20.12))") i, pebuf(i, 4)*to_mev_per_atom, pebuf(i, 1)*to_mev_per_atom, & + pebuf(i, 2)*to_mev_per_atom, pebuf(i, 3)*to_mev_per_atom + end do + + close (u) + + write (*, '(A)') ' ... energies writen to `outfile.energies`' + end if + end if + + +end block epotthings + +call mw%destroy() + +end program \ No newline at end of file diff --git a/src/effective_hamiltonian/manual.md b/src/effective_hamiltonian/manual.md new file mode 100644 index 00000000..e69de29b diff --git a/src/effective_hamiltonian/options.f90 b/src/effective_hamiltonian/options.f90 new file mode 100644 index 00000000..d8b52550 --- /dev/null +++ b/src/effective_hamiltonian/options.f90 @@ -0,0 +1,104 @@ +#include "precompilerdefinitions" +module options +use konstanter, only: flyt, r8, lo_hugeint, lo_huge, lo_author, lo_version, lo_licence, lo_tiny, lo_status, lo_exitcode_baddim, lo_exitcode_param +use gottochblandat, only: lo_stop_gracefully +use flap, only: command_line_interface +implicit none +private +public :: lo_opts + +type lo_opts + integer :: verbosity = -lo_hugeint + logical :: thirdorder = .false. + logical :: fourthorder = .false. + integer :: stride = 1 !! TODO CURRENTLY STRIDE IS IGNORED WHEN PARSING infile.positions + integer :: nconf = -1 + logical :: quantum = .false. + real(flyt) :: temperature = -lo_huge +contains + procedure :: parse +end type + +contains + +subroutine parse(opts) + !> the options + class(lo_opts), intent(out) :: opts + !> the helper parser + type(command_line_interface) :: cli + + logical :: dumlog + real(r8), dimension(3) :: dumr8v + integer :: errctr + + ! basic info + call cli%init(progname='effective_hamiltonian', & + authors=lo_author, & + version=lo_version, & + license=lo_licence, & + help='Usage: ', & + description='Calculates the TDEP effective Hamiltonian from force constants. & + &By default this command will attempt to load results form an MD simulation. & + &Either from the result of pack_simulation or the unpacked results of an MD simulation. & + &If the nconf flag is set then MD data will be ignored and configuratoins will be genereated & + &from a harmonic distribution according to the quantum and temperature flags.', & + examples=["mpirun effective_hamiltonian --thirdorder --fourthorder ", & + "mpirun effective_hamiltonian --thirdorder --fourthorder --nconf 200"], & + epilog=new_line('a')//"...") + + ! Specify some options + call cli%add(switch='--thirdorder', & + help='Compute third order contribution to the potential energy', & + required=.false., act='store_true', def='.false.', error=lo_status) + if (lo_status .ne. 0) stop + call cli%add(switch='--fourthorder', & + help='Compute fourth order contribution to the potential energy', & + required=.false., act='store_true', def='.false.', error=lo_status) + if (lo_status .ne. 0) stop + call cli%add(switch='--stride', switch_ab='s', & + help='Use every N configuration instead of all.', & + required=.false., act='store', def='1', error=lo_status) + if (lo_status .ne. 0) stop + call cli%add(switch='--nconf', switch_ab='-s', & + help='Automatically generate configurations, as done by the canonical configurations command. If not passed, an infile.positions is required.', & + required=.false., act='store', def='-1', error=lo_status) + if (lo_status .ne. 0) stop + call cli%add(switch='--quantum', & + help='If using --nconf, generate quantum or classical configurations.', & + required=.false., act='store_true', def='.false.', error=lo_status) + if (lo_status .ne. 0) stop + call cli%add(switch='--temperature', switch_ab='-t', & + help='If using --nconf, generate configurations at this temperature.', & + required=.false., act='store', def='-1', error=lo_status) + if (lo_status .ne. 0) stop + cli_manpage + cli_verbose + + ! actually parse it + errctr = 0 + call cli%parse(error=lo_status) + if (lo_status .ne. 0) stop + ! Should the manpage be generated? In that case, no point to go further than here. + dumlog = .false. + call cli%get(switch='--manpage', val=dumlog, error=lo_status); errctr = errctr + lo_status + if (dumlog) then + call cli%save_man_page(trim(cli%progname)//'.1') + call cli%save_usage_to_markdown(trim(cli%progname)//'.md') + write (*, *) 'Wrote manpage for "'//trim(cli%progname)//'"' + stop + end if + call cli%get(switch='--verbose', val=dumlog, error=lo_status); errctr = errctr + lo_status; if (dumlog) opts%verbosity = 2 + + ! get real options + call cli%get(switch='--thirdorder', val=opts%thirdorder, error=lo_status); errctr = errctr + lo_status + call cli%get(switch='--fourthorder', val=opts%fourthorder, error=lo_status); errctr = errctr + lo_status + call cli%get(switch="--stride", val = opts%stride, error = lo_status); errctr = errctr + lo_status + call cli%get(switch="--nconf", val = opts%nconf, error = lo_status); errctr = errctr + lo_status + call cli%get(switch="--temperature", val = opts%temperature, error = lo_status); errctr = errctr + lo_status + call cli%get(switch='--quantum', val = opts%quantum, error=lo_status); errctr = errctr + lo_status + + if (errctr .ne. 0) call lo_stop_gracefully(['Failed parsing the command line options'], lo_exitcode_baddim) + +end subroutine + +end module diff --git a/src/libolle/Makefile b/src/libolle/Makefile index fba9d326..22778fe2 100644 --- a/src/libolle/Makefile +++ b/src/libolle/Makefile @@ -140,7 +140,8 @@ $(OBJECT_PATH)type_forcemap_coefficient_diel.o \ $(OBJECT_PATH)type_forcemap_returntensors.o \ $(OBJECT_PATH)type_forcemap_io.o \ $(OBJECT_PATH)type_lassosolvers.o\ -$(OBJECT_PATH)lo_fftgrid_helper.o +$(OBJECT_PATH)lo_fftgrid_helper.o \ +$(OBJECT_PATH)lo_epot.o \ # maybe CGAL, maybe not. ifeq ($(USECGAL),yes) OBJScg = \ @@ -883,5 +884,19 @@ $(OBJECT_PATH)lo_fftgrid_helper.o:\ $(OBJECT_PATH)lo_randomnumbers.o $(FC) $(FFLAGS) -c lo_fftgrid_helper.f90 -o $@ +$(OBJECT_PATH)lo_epot.o:\ + lo_epot.f90 \ + $(OBJECT_PATH)konstanter.o\ + $(OBJECT_PATH)type_mdsim.o\ + $(OBJECT_PATH)gottochblandat.o\ + $(OBJECT_PATH)mpi_wrappers.o\ + $(OBJECT_PATH)lo_memtracker.o\ + $(OBJECT_PATH)type_crystalstructure.o\ + $(OBJECT_PATH)type_forceconstant_firstorder.o\ + $(OBJECT_PATH)type_forceconstant_secondorder.o\ + $(OBJECT_PATH)type_forceconstant_thirdorder.o\ + $(OBJECT_PATH)type_forceconstant_fourthorder.o + $(FC) $(FFLAGS) -c lo_epot.f90 -o $@ + clean: rm -rf $(OBJECT_PATH)*.o $(MODULE_PATH)*mod diff --git a/src/libolle/lo_epot.f90 b/src/libolle/lo_epot.f90 new file mode 100644 index 00000000..698fc2e0 --- /dev/null +++ b/src/libolle/lo_epot.f90 @@ -0,0 +1,503 @@ +module lo_epot + !! Deal with many kinds of potential energy differences + use konstanter, only: r8, lo_pi, lo_twopi, lo_tol, lo_sqtol, lo_status, lo_Hartree_to_eV, lo_kb_hartree + use gottochblandat, only: tochar, walltime, lo_chop, lo_trueNtimes, lo_progressbar_init, & + lo_progressbar, lo_frobnorm, open_file, lo_flattentensor, lo_sqnorm, lo_outerproduct, lo_mean, & + lo_points_on_sphere, lo_mean, lo_stddev + use mpi_wrappers, only: lo_mpi_helper, lo_stop_gracefully + use lo_memtracker, only: lo_mem_helper + use type_crystalstructure, only: lo_crystalstructure + use type_forceconstant_secondorder, only: lo_forceconstant_secondorder + use type_forceconstant_thirdorder, only: lo_forceconstant_thirdorder + use type_forceconstant_fourthorder, only: lo_forceconstant_fourthorder + use type_mdsim, only: lo_mdsim + implicit none + + private + public :: lo_energy_differences + + type lo_energy_differences + ! forceconstants needed for evaluation + type(lo_forceconstant_secondorder) :: fc2 + type(lo_forceconstant_thirdorder) :: fc3 + type(lo_forceconstant_fourthorder) :: fc4 + real(r8), dimension(:, :, :, :), allocatable :: fcp + logical :: polar = .false. + + contains + procedure :: setup => setup_potential_energy_differences + procedure :: statistical_sampling + procedure :: energies_and_forces + end type + + contains + + !> statistically sample + subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, quantum, ebuf, mw, mem, verbosity) + !> container for potential energy differences + class(lo_energy_differences), intent(inout) :: pot + !> unitcell + type(lo_crystalstructure), intent(in) :: uc + !> supercell + type(lo_crystalstructure), intent(in) :: ss + !> second order forceconstant + type(lo_forceconstant_secondorder), intent(in) :: fc + !> number of configs to generate + integer, intent(in) :: nstep + !> temperature + real(r8), intent(in) :: temperature + !> quantum configurations? + logical :: quantum + !> storage for potential energies + real(r8), dimension(:, :), intent(out) :: ebuf + !> MPI helper + type(lo_mpi_helper), intent(inout) :: mw + !> memory tracker + type(lo_mem_helper), intent(inout) :: mem + !> talk a lot? + integer, intent(in) :: verbosity + + type(lo_crystalstructure) :: p + integer :: ctr, i + real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp + real(r8) :: e2, e3, e4, ep + + ! Copy of structure to work with + p = ss + + ebuf = 0.0_r8 + + ! Dummy space for force + call mem%allocate(f2, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(f3, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(f4, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(fp, [3, ss%na], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + f2 = 0.0_r8 + f3 = 0.0_r8 + f4 = 0.0_r8 + fp = 0.0_r8 + + do i = 1, nstep + + if (mod(i, mw%n) .ne. mw%r) cycle + + ! Reset structure + p%u = 0.0_r8 + p%v = 0.0_r8 + p%r = ss%r + ! Get new structure + call pot%fc2%initialize_cell(p, uc, fc, temperature, quantum, .false., -1.0_r8, mw, nosync=.true.) + ! Calculate the energy + call pot%energies_and_forces(p%u, e2, e3, e4, ep, f2, f3, f4, fp) + + ! ek = p%kinetic_energy()/(p%na) + + ebuf(i, 1) = e2 + ebuf(i, 2) = e3 + ebuf(i, 3) = e4 + ebuf(i, 4) = ep + end do + + call mw%allreduce('sum', ebuf) + + call mem%deallocate(f2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%deallocate(f3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%deallocate(f4, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%deallocate(fp, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + end subroutine + + ! !> calculate energies and forces for a given configuration + subroutine energies_and_forces(pot, u, e2, e3, e4, ep, f2, f3, f4, fp) + !> container for potential energy differences + class(lo_energy_differences), intent(in) :: pot + !> displacements + real(r8), dimension(:, :), intent(in) :: u + !> energies + real(r8), intent(out) :: e2, e3, e4, ep + !> forces + real(r8), dimension(:, :), intent(out) :: f2, f3, f4, fp + + real(r8), dimension(3, 3, 3, 3) :: m4 + real(r8), dimension(3, 3, 3) :: m3 + real(r8), dimension(3, 3) :: m2 + real(r8), dimension(3) :: v0, u2, u3, u4 + integer :: a1, a2, a3, a4, i1, i2, i3, i4, i + + e2 = 0.0_r8 + e3 = 0.0_r8 + e4 = 0.0_r8 + ep = 0.0_r8 + f2 = 0.0_r8 + f3 = 0.0_r8 + f4 = 0.0_r8 + fp = 0.0_r8 + + ! pair term + do a1 = 1, size(u, 2) + v0 = 0.0_r8 + do i1 = 1, pot%fc2%atom(a1)%n + a2 = pot%fc2%atom(a1)%pair(i1)%i2 + m2 = pot%fc2%atom(a1)%pair(i1)%m + v0 = v0 - matmul(m2, u(:, a2)) + end do + e2 = e2 - dot_product(u(:, a1), v0)*0.5_r8 + f2(:, a1) = v0 + end do + + ! polar term + if (pot%polar) then + ep = 0.0_r8 + do a1 = 1, size(u, 2) + v0 = 0.0_r8 + do a2 = 1, size(u, 2) + v0 = v0 - matmul(pot%fcp(:, :, a1, a2), u(:, a2)) + end do + ep = ep - dot_product(u(:, a1), v0)*0.5_r8 + fp(:, a1) = v0 + end do + end if + + ! triplet term + if (pot%fc3%na .eq. size(u, 2)) then + do a1 = 1, size(u, 2) + v0 = 0.0_r8 + do i = 1, pot%fc3%atom(a1)%n + m3 = pot%fc3%atom(a1)%triplet(i)%m + a2 = pot%fc3%atom(a1)%triplet(i)%i2 + a3 = pot%fc3%atom(a1)%triplet(i)%i3 + u2 = u(:, a2) + u3 = u(:, a3) + do i1 = 1, 3 + do i2 = 1, 3 + do i3 = 1, 3 + v0(i1) = v0(i1) - m3(i1, i2, i3)*u2(i2)*u3(i3) + end do + end do + end do + end do + v0 = v0*0.5_r8 + f3(:, a1) = v0 + e3 = e3 - dot_product(v0, u(:, a1))/3.0_r8 + end do + end if + + ! quartet term + if (pot%fc4%na .eq. size(u, 2)) then + do a1 = 1, size(u, 2) + v0 = 0.0_r8 + do i = 1, pot%fc4%atom(a1)%n + m4 = pot%fc4%atom(a1)%quartet(i)%m + a2 = pot%fc4%atom(a1)%quartet(i)%i2 + a3 = pot%fc4%atom(a1)%quartet(i)%i3 + a4 = pot%fc4%atom(a1)%quartet(i)%i4 + u2 = u(:, a2) + u3 = u(:, a3) + u4 = u(:, a4) + do i1 = 1, 3 + do i2 = 1, 3 + do i3 = 1, 3 + do i4 = 1, 3 + v0(i1) = v0(i1) - m4(i1, i2, i3, i4)*u2(i2)*u3(i3)*u4(i4) + end do + end do + end do + end do + end do + v0 = v0/6.0_r8 + f4(:, a1) = v0 + e4 = e4 - dot_product(v0, u(:, a1))/4.0_r8 + end do + end if + end subroutine + + !> Calculate potential energy differences in several ways + subroutine setup_potential_energy_differences(pot, uc, ss, fc2, fc3, fc4, mw, verbosity) + !> container for potential energy differences + class(lo_energy_differences), intent(out) :: pot + !> unitcell + type(lo_crystalstructure), intent(inout) :: uc + !> supercell + type(lo_crystalstructure), intent(inout) :: ss + !> second order forceconstant + type(lo_forceconstant_secondorder), intent(in) :: fc2 + !> third order forceconstant + type(lo_forceconstant_thirdorder), intent(in) :: fc3 + !> fourth order forceconstant + type(lo_forceconstant_fourthorder), intent(in) :: fc4 + !> mpi things + type(lo_mpi_helper), intent(inout) :: mw + !> how much to talk + integer, intent(in) :: verbosity + + real(r8) :: timer, t0, t1 + + timer = walltime() + t0 = timer + t1 = timer + + if (verbosity .gt. 0) then + write (*, *) '' + write (*, *) 'PREPARING POTENTIAL ENERGY DIFFERENCES' + end if + + call fc2%remap(uc, ss, pot%fc2) + if (verbosity .gt. 0) then + t1 = walltime() + write (*, *) '... remapped second order (', tochar(t1 - t0), ')' + t0 = t1 + end if + + if (fc3%na .gt. 0) then + call fc3%remap(uc, ss, pot%fc3) + if (verbosity .gt. 0) then + t1 = walltime() + write (*, *) '... remapped third order (', tochar(t1 - t0), ')' + t0 = t1 + end if + end if + + if (fc4%na .gt. 0) then + call fc4%remap(uc, ss, pot%fc4) + if (verbosity .gt. 0) then + t1 = walltime() + write (*, *) '... remapped fourth order (', tochar(t1 - t0), ')' + t0 = t1 + end if + end if + + if (fc2%polar) then + allocate (pot%fcp(3, 3, ss%na, ss%na)) + pot%fcp = 0.0_r8 + call fc2%supercell_longrange_dynamical_matrix_at_gamma(ss, pot%fcp, 1E-15_r8) + pot%polar = .true. + else + pot%polar = .false. + end if + if (verbosity .gt. 0) then + t1 = walltime() + write (*, *) '... built polar forceconstant (', tochar(t1 - t0), ')' + t0 = t1 + end if + end subroutine + + ! !> remove longrange interactions + ! subroutine potential_energy_difference(sim,uc,ss,fc,fct,fcf,Uref,mw,verbosity) + ! !> force-displacement data + ! type(lo_mdsim), intent(inout) :: sim + ! !> unitcell + ! type(lo_crystalstructure), intent(inout) :: uc + ! !> supercell + ! type(lo_crystalstructure), intent(inout) :: ss + ! !> second order forceconstant + ! type(lo_forceconstant_secondorder), intent(in) :: fc + ! !> third order forceconstant + ! type(lo_forceconstant_thirdorder), intent(in) :: fct + ! !> fourth order forceconstant + ! type(lo_forceconstant_fourthorder), intent(in) :: fcf + ! !> reference energy + ! real(flyt), intent(in) :: Uref + ! !> mpi things + ! type(lo_mpi_helper), intent(inout) :: mw + ! !> how much to talk + ! integer, intent(in) :: verbosity + ! + ! type(lo_forceconstant_secondorder) :: fcss + ! type(lo_forceconstant_thirdorder) :: fctss + ! type(lo_forceconstant_fourthorder) :: fcfss + ! real(flyt), dimension(:,:,:,:), allocatable :: forceconstant + ! real(flyt), dimension(:,:), allocatable :: fp,f2,f3,f4,de + ! ! real(flyt), dimension(:,:), allocatable :: f,de + ! ! real(flyt), dimension(3,3) :: m0,m1 + ! ! real(flyt), dimension(3) :: v0 + ! ! real(flyt) :: epot,epot_2,epot_3,epot_4,epot_dd + ! ! integer :: a1,a2,t,u + ! + ! init: block + ! real(flyt), dimension(3,3) :: m0,m1 + ! real(flyt) :: t0,f0 + ! integer :: a1,a2 + ! if ( mw%talk ) then + ! write(*,*) '' + ! write(*,*) 'CALCULATING POTENTIAL ENERGY DIFFERENCES' + ! endif + ! ! Remap the forceconstants to the supercell + ! call ss%classify('supercell',uc) + ! call fc%remap(uc,ss,fcss) + ! call fct%remap(uc,ss,fctss) + ! call fcf%remap(uc,ss,fcfss) + ! if ( mw%talk ) write(*,*) '... remapped forceconstants' + ! ! Get the dipole-dipole forceconstant + ! if ( fc%polar ) then + ! lo_allocate(forceconstant(3,3,ss%na,ss%na)) + ! call fc%supercell_longrange_dynamical_matrix_at_gamma(ss,forceconstant,1E-15_flyt) + ! ! small sanity-check, just in case + ! f0=0.0_flyt + ! do a1=1,ss%na + ! do a2=a1,ss%na + ! m0=forceconstant(:,:,a1,a2) + ! m1=transpose(forceconstant(:,:,a2,a1)) + ! f0=f0+lo_frobnorm(m0-m1) + ! enddo + ! enddo + ! f0=f0/(ss%na**2) + ! if ( f0 .gt. lo_tol ) then + ! write(*,*) 'ERROR, non-Hermitian electrostatic forceconstant: ',f0 + ! write(*,*) 'If this keeps happening I should not hard-code the Ewald tolerance.' + ! stop + ! endif + ! endif + ! + ! ! some space for energies and forces + ! lo_allocate(de(sim%nt,5)) + ! lo_allocate(fp(3,ss%na)) + ! lo_allocate(f2(3,ss%na)) + ! lo_allocate(f3(3,ss%na)) + ! lo_allocate(f4(3,ss%na)) + ! de=0.0_flyt + ! fp=0.0_flyt + ! f2=0.0_flyt + ! f3=0.0_flyt + ! f4=0.0_flyt + ! end block init + ! + ! forceenergy: block + ! real(flyt), dimension(3,3,3,3) :: m4 + ! real(flyt), dimension(3,3,3) :: m3 + ! real(flyt), dimension(3,3) :: m2 + ! real(flyt), dimension(3) :: u1,u2,u3,u4,v0 + ! real(flyt) :: energy + ! integer :: a1,a2,a3,a4,i1,i2,i3,i4,i,j,k,l,t + ! + ! do t=1,sim%nt + ! ! get the raw energy + ! de(t,1)=sim%stat%potential_energy(t) + ! + ! ! secondorder energy, first the polar term (if applicable) + ! energy=0.0_flyt + ! if ( fc%polar ) then + ! do a1=1,ss%na + ! v0=0.0_flyt + ! do a2=1,ss%na + ! v0=v0+matmul(forceconstant(:,:,a2,a1),u(:,a2,t)) + ! enddo + ! energy=energy+dot_product(u(:,a1,t),v0)*0.5_flyt + ! enddo + ! endif + ! ! then the pair term, always there + ! do a1=1,ss%na + ! v0=0.0_flyt + ! do i1=1,fcss%atom(a1)%n + ! a2=fcss%atom(a1)%pair(i1)%i2 + ! m2=fcss%atom(a1)%pair(i1)%m + ! v0=v0+matmul(m2,u(:,a2,t)) + ! enddo + ! energy=energy+dot_product(u(:,a1,t),v0)*0.5_flyt + ! enddo + ! de(t,2)=energy + ! + ! ! triplet term + ! f3=0.0_flyt + ! energy=0.0_flyt + ! do a1=1,fcss%na + ! v0=0.0_flyt + ! do i=1,fctss%atom(a1)%n + ! m3=fctss%atom(a1)%triplet(i)%m + ! a2=fctss%atom(a1)%triplet(i)%i2 + ! a3=fctss%atom(a1)%triplet(i)%i3 + ! u2=u(:,a2,t) + ! u3=u(:,a3,t) + ! do i1=1,3 + ! do i2=1,3 + ! do i3=1,3 + ! v0(i1)=v0(i1)-m3(i1,i2,i3)*u2(i2)*u3(i3) + ! enddo + ! enddo + ! enddo + ! enddo + ! f3(:,a1)=v0*0.5_flyt + ! energy=energy-dot_product(f3(:,a1),u(:,a1,t))/3.0_flyt + ! enddo + ! de(t,3)=energy + ! + ! ! quartet term + ! f4=0.0_flyt + ! energy=0.0_flyt + ! do a1=1,fcfss%na + ! v0=0.0_flyt + ! do i=1,fcfss%atom(a1)%n + ! m4=fcfss%atom(a1)%quartet(i)%m + ! a2=fcfss%atom(a1)%quartet(i)%i2 + ! a3=fcfss%atom(a1)%quartet(i)%i3 + ! a4=fcfss%atom(a1)%quartet(i)%i4 + ! u2=u(:,a2,t) + ! u3=u(:,a3,t) + ! u4=u(:,a4,t) + ! do i1=1,3 + ! do i2=1,3 + ! do i3=1,3 + ! do i4=1,3 + ! v0(i1)=v0(i1)-m4(i1,i2,i3,i4)*u2(i2)*u3(i3)*u4(i4) + ! enddo + ! enddo + ! enddo + ! enddo + ! enddo + ! f4(:,a1)=v0/6.0_flyt + ! energy=energy-dot_product(f4(:,a1),u(:,a1,t))/4.0_flyt + ! enddo + ! de(t,4)=energy + ! enddo + ! end block forceenergy + ! + ! if ( mw%talk ) then + ! write(*,*) ' raw avg potential energy:',lo_mean(de(:,1)) + ! write(*,*) ' E-E2:',lo_mean(de(:,1)-de(:,2)) + ! write(*,*) ' E-E2-E3:',lo_mean(de(:,1)-de(:,2)-de(:,3)) + ! write(*,*) ' E-E2-E3-E4:',lo_mean(de(:,1)-de(:,2)-de(:,3)-de(:,4)) + ! endif + ! + ! + ! !! Energy per atom + ! !de=1000*de/sim%na + ! ! + ! !if ( mw%talk ) then + ! !write(*,*) 'RAW:',lo_mean(de(:,1)),lo_stddev(de(:,1)) + ! !write(*,*) 'RAW:',lo_mean(de(:,1)-de(:,3)),lo_stddev(de(:,1)-de(:,3)) + ! !write(*,*) 'RAW:',lo_mean(de(:,1)-de(:,2)-de(:,3)),lo_stddev(de(:,1)-de(:,2)-de(:,3)) + ! !write(*,*) 'RAW:',lo_mean(de(:,1)-de(:,2)-de(:,3)-de(:,4)),lo_stddev(de(:,1)-de(:,2)-de(:,3)-de(:,4)) + ! !write(*,*) 'RAW:',lo_mean(de(:,1)-de(:,2)-de(:,3)-de(:,4)-de(:,5)),lo_stddev(de(:,1)-de(:,2)-de(:,3)-de(:,4)-de(:,5)) + ! !u=open_file('out','outfile.U0') + ! ! do t=1,sim%nt + ! ! write(u,"(6(1X,E18.12))") de(t,:)/1000*sim%na,(de(t,1)-de(t,2)-de(t,3)-de(t,4)-de(t,5))/1000*sim%na + ! ! enddo + ! !close(u) + ! !endif + ! + ! ! t0=walltime() + ! ! call lo_progressbar_init() + ! ! ! Subtract forces + ! ! lo_allocate(f(3,ss%na)) + ! ! do t=1,sim%nt + ! ! ! calculate the long-range forces + ! ! f=0.0_flyt + ! ! do a1=1,ss%na + ! ! do a2=1,ss%na + ! ! v0=matmul(forceconstant(:,:,a1,a2),u(:,a2,t)) + ! ! f(:,a1)=f(:,a1)-v0 + ! ! enddo + ! ! enddo + ! ! ! sanity check that they add up to zero + ! ! if ( abs(sum(f)) .gt. lo_sqtol ) then + ! ! write(*,*) '' + ! ! write(*,*) 'ERROR: dipole-dipole forces do not add up to zero!' + ! ! stop + ! ! endif + ! ! ! subtract them + ! ! sim%f(:,:,t)=sim%f(:,:,t)-f + ! ! if( lo_trueNtimes(t,50,sim%nt) ) call lo_progressbar(' ... subtracting dipole-dipole forces',t,sim%nt) + ! ! enddo + ! ! call lo_progressbar(' ... subtracting dipole-dipole forces',sim%nt,sim%nt,walltime()-t0) + ! end subroutine + + end module + \ No newline at end of file