From 105e444b9dd65319b6a82c7ef780c073ae87d5ae Mon Sep 17 00:00:00 2001 From: ejmeitz Date: Fri, 28 Mar 2025 11:01:53 -0400 Subject: [PATCH 1/8] move epot to libolle, start effective_hamiltonian Conflicts: src/anharmonic_free_energy/main.f90 --- src/anharmonic_free_energy/Makefile | 5 +- src/anharmonic_free_energy/epot.f90 | 534 -------------------------- src/anharmonic_free_energy/main.f90 | 2 +- src/effective_hamiltonian/Makefile | 31 ++ src/effective_hamiltonian/main.f90 | 104 +++++ src/effective_hamiltonian/manual.md | 0 src/effective_hamiltonian/options.f90 | 99 +++++ src/libolle/Makefile | 17 +- src/libolle/lo_epot.f90 | 508 ++++++++++++++++++++++++ 9 files changed, 761 insertions(+), 539 deletions(-) delete mode 100644 src/anharmonic_free_energy/epot.f90 create mode 100644 src/effective_hamiltonian/Makefile create mode 100644 src/effective_hamiltonian/main.f90 create mode 100644 src/effective_hamiltonian/manual.md create mode 100644 src/effective_hamiltonian/options.f90 create mode 100644 src/libolle/lo_epot.f90 diff --git a/src/anharmonic_free_energy/Makefile b/src/anharmonic_free_energy/Makefile index 2eeeeee8..2d288a40 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) @@ -30,6 +30,5 @@ $(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..f085e247 --- /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 = -lolle -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 $(LIBS) -o $@ +$(OBJECT_PATH)options.o: + $(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@ + diff --git a/src/effective_hamiltonian/main.f90 b/src/effective_hamiltonian/main.f90 new file mode 100644 index 00000000..acdf6955 --- /dev/null +++ b/src/effective_hamiltonian/main.f90 @@ -0,0 +1,104 @@ +#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 +use mpi_wrappers, only: lo_mpi_helper +use lo_memtracker, only: lo_mem_helper +use gottochblandat, only: tochar, walltime, lo_stop_gracefully, open_file + +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: statistical_sampling, setup_potential_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 + +logical :: generate_configs = .false. + +init: block + ! Get CLI options + call opts%parse() + + 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 scattering ', opts%thirdorder + write (*, '(1X,A40,L3)') 'Fourthorder scattering ', opts%fourthorder + write (*, '(1X,A40,I5)') 'Stride ', opts%stride + write(*, '(1X,A40,L3)') 'Generate canonical configs ', generate_configs + write (*, '(1X,A40,L3)') 'Quantum configurations ', opts%quantum + write (*, '(1X,A40,F20.12)') 'Temperature ', opts%temperature + end if + + + call mw%init() + call mem%init() + if (.not. mw%talk) opts%verbosity = -100 + + ! 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 fc%readfromfile(uc, 'infile.forceconstant', mem, verbosity=-1) + if (opts%thirdorder) call fct%readfromfile(uc, 'infile.forceconstant_thirdorder') + if (opts%fourthorder) call fcf%readfromfile(uc, 'infile.forceconstant_fourthorder') + if (mw%talk) write (*, *) '... read forceconstants' + + if (generate_configs) then + if (opts%temperature .lt. 0.0_r8) then + call lo_stop_gracefully(['You need to provide a positive temperature to generate configurations'], & + lo_exitcode_param, __FILE__, __LINE__, mw%comm) + end if + if (mw%talk .and. opts%quantum) write (*, *) '... will generate quantum configurations' + if (mw%talk .and. opts%quantum) write (*, *) '... will generate classical configurations' + end if + + call pot%setup(uc, ss, fc2, fc3, fc4, mw, verbosity) + if (mw%talk) write (*, *) '... setup potential energy calculator' + + if (.not. generate_configs) then + readmag = .false. + readdiel = .false. + call sim%read_from_file(verbosity=opts%verbosity + 2, stride=opts%stride, mw=mw) + if (mw%talk) write (*, *) '... parsed infile.positions' + end if + +end block init + + +energy : block + + real(r8), dimension(:, :), allocatable :: ebuf + + if (generate_configs) then + call pot%statistical_sampling(uc, ss, fc2, otps%nconf, opts%temperature, mw, mem, opts%verbosity) + else + ! use parsed sim to evalulate energies + end if + + +end block energy \ 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..9253ce8c --- /dev/null +++ b/src/effective_hamiltonian/options.f90 @@ -0,0 +1,99 @@ +#include "precompilerdefinitions" +module options +use konstanter, only: 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 + 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 only the potential energy component is calculated.', & + examples=["mpirun effective_hamiltonian --thirdorder --fourthorder"], & + epilog=new_line('a')//"...") + + ! Specify some options + call cli%add(switch='--thirdorder', & + help='Compute third order anharmonic correction to the free 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 anharmonic correction to the free 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..605c8780 --- /dev/null +++ b/src/libolle/lo_epot.f90 @@ -0,0 +1,508 @@ +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, nstep, temperature, 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 + !> 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 :: nstep, ctr, i + real(r8), dimension(:, :), allocatable, intent(out) :: ebuf + real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp + real(r8) :: e2, e3, e4, ep, 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, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + 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 + + tomev = lo_Hartree_to_eV*1000 + + ctr = 0 + do i = 1, nstep + + ctr = ctr + 1 ! does parallelizing like this create a copy of the IFCs on each rank? not ideal + if (mod(ctr, 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, .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) + + ebuf(i, 1) = e2 + ebuf(i, 2) = e3 + ebuf(i, 3) = e4 + ebuf(i, 4) = ep + end do + + call mw%barrier() ! needed? + + 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 From b17851ff92bfe576ee776d36d0752ad48613a93d Mon Sep 17 00:00:00 2001 From: ejmeitz Date: Fri, 28 Mar 2025 14:12:24 -0400 Subject: [PATCH 2/8] main code builds --- build_things.sh | 1 + src/effective_hamiltonian/main.f90 | 171 +++++++++++++++++++++++--- src/effective_hamiltonian/options.f90 | 10 +- src/libolle/lo_epot.f90 | 22 ++-- 4 files changed, 164 insertions(+), 40 deletions(-) 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/effective_hamiltonian/main.f90 b/src/effective_hamiltonian/main.f90 index acdf6955..cdf22d78 100644 --- a/src/effective_hamiltonian/main.f90 +++ b/src/effective_hamiltonian/main.f90 @@ -2,17 +2,18 @@ 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_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 +use gottochblandat, only: tochar, walltime, lo_stop_gracefully, open_file, lo_progressbar_init, & + lo_progressbar, lo_does_file_exist, lo_trueNtimes 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: statistical_sampling, setup_potential_energy_differences +use lo_epot, only: lo_energy_differences use type_crystalstructure, only: lo_crystalstructure use options, only: lo_opts @@ -30,10 +31,16 @@ program effective_hamiltonian type(lo_mpi_helper) :: mw type(lo_mem_helper) :: mem +real(r8), dimension(:, :), allocatable :: ebuf logical :: generate_configs = .false. init: block + + integer :: f, i, j, l, readrank + logical :: readonthisrank, mpiparallel + real(r8) :: t0 + ! Get CLI options call opts%parse() @@ -41,8 +48,8 @@ program effective_hamiltonian if (mw%talk) then write (*, *) 'Recap of the parameters governing the calculation' - write (*, '(1X,A40,L3)') 'Thirdorder scattering ', opts%thirdorder - write (*, '(1X,A40,L3)') 'Fourthorder scattering ', opts%fourthorder + 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 write (*, '(1X,A40,L3)') 'Quantum configurations ', opts%quantum @@ -63,42 +70,166 @@ program effective_hamiltonian call uc%classify('spacegroup', timereversal=.true.) call ss%classify('supercell', uc) - call fc%readfromfile(uc, 'infile.forceconstant', mem, verbosity=-1) - if (opts%thirdorder) call fct%readfromfile(uc, 'infile.forceconstant_thirdorder') - if (opts%fourthorder) call fcf%readfromfile(uc, 'infile.forceconstant_fourthorder') + 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%temperature .lt. 0.0_r8) then - call lo_stop_gracefully(['You need to provide a positive temperature to generate configurations'], & - lo_exitcode_param, __FILE__, __LINE__, mw%comm) + 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' end if - call pot%setup(uc, ss, fc2, fc3, fc4, mw, verbosity) + call pot%setup(uc, ss, fc2, fc3, fc4, mw, opts%verbosity) if (mw%talk) write (*, *) '... setup potential energy calculator' if (.not. generate_configs) then - readmag = .false. - readdiel = .false. - call sim%read_from_file(verbosity=opts%verbosity + 2, stride=opts%stride, mw=mw) - if (mw%talk) write (*, *) '... parsed infile.positions' + ! 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 + ! Just read the infile.positions/ infile.meta dont need forces or stat files + ! Will need to populate some fields of sim manually + + ! should I just do sim%readfromfile? would let me calculate things like U0 and other cumulants + ! to get infile.positions you need to run MD anyway and might as well dump the others + + sim%crystalstructure = ss ! pretty sure this copies, but oh well its not a lot of data + + + ! We are reading on one rank, then broadcasting + readrank = 0 + if (mw%r .eq. readrank) then + readonthisrank = .true. + else + readonthisrank = .false. + end if + + + if (readonthisrank) then + f = open_file('in', 'infile.meta') + read (f, *) sim%na + read (f, *) sim%nt + close (f) + end if + + lo_allocate(sim%r(3, sim%na, sim%nt)) + sim%r = 0.0_r8 + + if (readonthisrank) then + t0 = walltime() + if (opts%verbosity .gt. 0) call lo_progressbar_init() + f = open_file('in', 'infile.positions') + i = 0 + do l = 1, sim%nt + + i = i + 1 + do j = 1, sim%na + read (f, *) sim%r(:, j, i) + end do + + if (opts%verbosity .gt. 0) then + if (lo_trueNtimes(l, 25, sim%nt)) then + call lo_progressbar('... reading positions', l, sim%nt, walltime() - t0) + end if + end if + end do + close (f) + end if + + call mw%bcast(sim%r, readrank, __FILE__, __LINE__) + call mw%bcast(sim%na, readrank, __FILE__, __LINE__) + call mw%bcast(sim%nt, readrank, __FILE__, __LINE__) + + ! get non-pbc positions + call sim%get_nonpbc_positions() + if (opts%verbosity .gt. 0) write (*, *) '... unwrapped positions' + ! get displacements + call sim%get_displacements() + if (opts%verbosity .gt. 0) write (*, *) '... got displacements' + + end if + if (mw%talk) write (*, *) '... parsed simulation data' end if end block init + + energy : block - real(r8), dimension(:, :), allocatable :: ebuf + integer :: ctr, i + real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp + real(r8) :: e2, e3, e4, ep + + call mem%allocate(ebuf, [opts%nconf, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + ebuf = 0.0_r8 if (generate_configs) then - call pot%statistical_sampling(uc, ss, fc2, otps%nconf, opts%temperature, mw, mem, opts%verbosity) + call pot%statistical_sampling(uc, ss, fc2, opts%nconf, opts%temperature, ebuf, mw, mem, opts%verbosity) else - ! use parsed sim to evalulate energies + + ! 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? + ctr = 0 + do i = 1, sim%nt !! TODO SKIP BASED ON STRIDE?? + + ctr = ctr + 1 + if (mod(ctr, 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) + + ebuf(i, 1) = e2 + ebuf(i, 2) = e3 + ebuf(i, 3) = e4 + ebuf(i, 4) = ep + end do end if + if (mw%talk) write (*, *) '... calculated energies' + + +end block energy + + +io : block + + integer :: i, u + real(r8) :: total_energy + + ! file for the energies + u = open_file('out', 'outfile.energies') + write (u, '(A,A)') '# Unit: ', 'eV/atom' + write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) + write (u, "(A)") '# conf Etotal Epolar & + &Epair Etriplet Equartet' + + do i = 1, sim%nt + total_energy = sum(ebuf(i,:)) + write (u, "(1X,I5,5(2X,E20.12))") i, total_energy*lo_Hartree_to_eV, ebuf(i, 4)*lo_Hartree_to_eV, & + ebuf(i, 1)*lo_Hartree_to_eV, ebuf(i, 2)*lo_Hartree_to_eV, & + ebuf(i, 3)*lo_Hartree_to_eV + end do + + ! close outfile.energies + write (*, '(A)') ' ... energies writen to `outfile.energies`' + close (u) + +end block io +call mw%destroy() -end block energy \ No newline at end of file +end program \ No newline at end of file diff --git a/src/effective_hamiltonian/options.f90 b/src/effective_hamiltonian/options.f90 index 9253ce8c..93a196dd 100644 --- a/src/effective_hamiltonian/options.f90 +++ b/src/effective_hamiltonian/options.f90 @@ -1,6 +1,6 @@ #include "precompilerdefinitions" module options -use konstanter, only: r8, lo_hugeint, lo_huge, lo_author, lo_version, lo_licence, lo_tiny, lo_status, lo_exitcode_baddim, lo_exitcode_param +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 @@ -11,7 +11,7 @@ module options integer :: verbosity = -lo_hugeint logical :: thirdorder = .false. logical :: fourthorder = .false. - integer :: stride = 1 + integer :: stride = 1 !! TODO CURRENTLY STRIDE IS IGNORED WHEN PARSING infile.positions integer :: nconf = -1 logical :: quantum = .false. real(flyt) :: temperature = -lo_huge @@ -50,11 +50,11 @@ subroutine parse(opts) help='Compute fourth order anharmonic correction to the free 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' & + 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' & + 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 @@ -62,7 +62,7 @@ subroutine parse(opts) 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' & + 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 diff --git a/src/libolle/lo_epot.f90 b/src/libolle/lo_epot.f90 index 605c8780..ef5f154e 100644 --- a/src/libolle/lo_epot.f90 +++ b/src/libolle/lo_epot.f90 @@ -1,4 +1,4 @@ -module epot +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, & @@ -6,7 +6,6 @@ module epot 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 @@ -55,7 +54,7 @@ subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, m integer, intent(in) :: verbosity type(lo_crystalstructure) :: p - integer :: nstep, ctr, i + integer :: ctr, i real(r8), dimension(:, :), allocatable, intent(out) :: ebuf real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp real(r8) :: e2, e3, e4, ep, tomev, emean @@ -63,9 +62,6 @@ subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, m ! 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, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) ebuf = 0.0_r8 @@ -79,9 +75,7 @@ subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, m f3 = 0.0_r8 f4 = 0.0_r8 fp = 0.0_r8 - - tomev = lo_Hartree_to_eV*1000 - + ctr = 0 do i = 1, nstep @@ -97,13 +91,11 @@ subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, m ! Calculate the energy call pot%energies_and_forces(p%u, e2, e3, e4, ep, f2, f3, f4, fp) - ebuf(i, 1) = e2 - ebuf(i, 2) = e3 - ebuf(i, 3) = e4 - ebuf(i, 4) = ep + ebuf(i, 1) = e2 / ss%na + ebuf(i, 2) = e3 / ss%na + ebuf(i, 3) = e4 / ss%na + ebuf(i, 4) = ep / ss%na end do - - call mw%barrier() ! needed? call mem%deallocate(f2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) call mem%deallocate(f3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) From 1da38b5f856fcea552b11b020c4cff2427e8e574 Mon Sep 17 00:00:00 2001 From: ejmeitz Date: Fri, 28 Mar 2025 18:06:12 -0400 Subject: [PATCH 3/8] it works! --- src/effective_hamiltonian/main.f90 | 215 ++++++++++++++------------ src/effective_hamiltonian/options.f90 | 13 +- src/libolle/lo_epot.f90 | 8 +- 3 files changed, 133 insertions(+), 103 deletions(-) diff --git a/src/effective_hamiltonian/main.f90 b/src/effective_hamiltonian/main.f90 index cdf22d78..7cb9e3a4 100644 --- a/src/effective_hamiltonian/main.f90 +++ b/src/effective_hamiltonian/main.f90 @@ -6,7 +6,7 @@ program effective_hamiltonian 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_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 @@ -31,7 +31,8 @@ program effective_hamiltonian type(lo_mpi_helper) :: mw type(lo_mem_helper) :: mem -real(r8), dimension(:, :), allocatable :: ebuf +real(r8), dimension(:, :), allocatable :: pebuf +! real(r8), dimension(:), allocatable :: kebuf logical :: generate_configs = .false. @@ -83,7 +84,7 @@ program effective_hamiltonian if (mw%talk .and. opts%quantum) write (*, *) '... will generate classical configurations' end if - call pot%setup(uc, ss, fc2, fc3, fc4, mw, opts%verbosity) + 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 @@ -91,66 +92,8 @@ program effective_hamiltonian 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 - ! Just read the infile.positions/ infile.meta dont need forces or stat files - ! Will need to populate some fields of sim manually - - ! should I just do sim%readfromfile? would let me calculate things like U0 and other cumulants - ! to get infile.positions you need to run MD anyway and might as well dump the others - - sim%crystalstructure = ss ! pretty sure this copies, but oh well its not a lot of data - - - ! We are reading on one rank, then broadcasting - readrank = 0 - if (mw%r .eq. readrank) then - readonthisrank = .true. - else - readonthisrank = .false. - end if - - - if (readonthisrank) then - f = open_file('in', 'infile.meta') - read (f, *) sim%na - read (f, *) sim%nt - close (f) - end if - - lo_allocate(sim%r(3, sim%na, sim%nt)) - sim%r = 0.0_r8 - - if (readonthisrank) then - t0 = walltime() - if (opts%verbosity .gt. 0) call lo_progressbar_init() - f = open_file('in', 'infile.positions') - i = 0 - do l = 1, sim%nt - - i = i + 1 - do j = 1, sim%na - read (f, *) sim%r(:, j, i) - end do - - if (opts%verbosity .gt. 0) then - if (lo_trueNtimes(l, 25, sim%nt)) then - call lo_progressbar('... reading positions', l, sim%nt, walltime() - t0) - end if - end if - end do - close (f) - end if - - call mw%bcast(sim%r, readrank, __FILE__, __LINE__) - call mw%bcast(sim%na, readrank, __FILE__, __LINE__) - call mw%bcast(sim%nt, readrank, __FILE__, __LINE__) - - ! get non-pbc positions - call sim%get_nonpbc_positions() - if (opts%verbosity .gt. 0) write (*, *) '... unwrapped positions' - ! get displacements - call sim%get_displacements() - if (opts%verbosity .gt. 0) write (*, *) '... got displacements' - + 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 @@ -162,17 +105,26 @@ program effective_hamiltonian energy : block - integer :: ctr, i + integer :: i, u real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp - real(r8) :: e2, e3, e4, ep + real(r8) :: e2, e3, e4, ep, total_energy, to_ev_per_atom - call mem%allocate(ebuf, [opts%nconf, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) - ebuf = 0.0_r8 if (generate_configs) then - call pot%statistical_sampling(uc, ss, fc2, opts%nconf, opts%temperature, ebuf, mw, mem, opts%verbosity) + + call mem%allocate(pebuf, [opts%nconf, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + ! call mem%allocate(kebuf, [opts%nconf], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + pebuf = 0.0_r8 + ! kebuf = 0.0_r8 + + call pot%statistical_sampling(uc, ss, fc2, opts%nconf, opts%temperature, pebuf, mw, mem, opts%verbosity) else + call mem%allocate(pebuf, [sim%nt, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + ! call mem%allocate(kebuf, [sim%nt], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + pebuf = 0.0_r8 + ! 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__) @@ -184,51 +136,124 @@ program effective_hamiltonian fp = 0.0_r8 ! how do I add a progress bar here without race condition? - ctr = 0 - do i = 1, sim%nt !! TODO SKIP BASED ON STRIDE?? + do i = 1, sim%nt - ctr = ctr + 1 - if (mod(ctr, mw%n) .ne. mw%r) cycle + 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) - ebuf(i, 1) = e2 - ebuf(i, 2) = e3 - ebuf(i, 3) = e4 - ebuf(i, 4) = ep + pebuf(i, 1) = e2 + pebuf(i, 2) = e3 + pebuf(i, 3) = e4 + pebuf(i, 4) = ep end do end if - if (mw%talk) write (*, *) '... calculated energies' + call mw%allreduce('sum', pebuf) -end block energy - + if (mw%talk) write (*, *) '... calculated energies' -io : block +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, to_ev_per_atom integer :: i, u - real(r8) :: total_energy - ! file for the energies - u = open_file('out', 'outfile.energies') - write (u, '(A,A)') '# Unit: ', 'eV/atom' - write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) - write (u, "(A)") '# conf Etotal Epolar & - &Epair Etriplet Equartet' + if (generate_configs) then + T_actual = opts%temperature + else + T_actual = sim%temperature_thermostat + end if + + ! Calculate the baseline energy + allocate (ediff(sim%nt, 5)) + ediff = 0.0_r8 do i = 1, sim%nt - total_energy = sum(ebuf(i,:)) - write (u, "(1X,I5,5(2X,E20.12))") i, total_energy*lo_Hartree_to_eV, ebuf(i, 4)*lo_Hartree_to_eV, & - ebuf(i, 1)*lo_Hartree_to_eV, ebuf(i, 2)*lo_Hartree_to_eV, & - ebuf(i, 3)*lo_Hartree_to_eV + 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 - ! close outfile.energies - write (*, '(A)') ' ... energies writen to `outfile.energies`' - close (u) + ! And normalize it to be per atom + cumulant = cumulant/real(ss%na, r8) + 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)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 1)*lo_Hartree_to_eV + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2): ', cumulant(1, 2)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 2)*lo_Hartree_to_eV + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar): ', cumulant(1, 3)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 3)*lo_Hartree_to_eV + if (opts%thirdorder) then + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar-E3): ', cumulant(1, 4)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 4)*lo_Hartree_to_eV + 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)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 5)*lo_Hartree_to_eV + 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)*lo_Hartree_to_eV + else if(opts%fourthorder) then + U0 = cumulant(1,5)*lo_Hartree_to_eV + else + U0 = cumulant(1,3)*lo_Hartree_to_eV + 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: ', 'eV/atom' + write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) + write (u, *) 'U0 [eV/atom]: ', U0 + write (u, "(A)") '# conf Etotal_actual Etotal_tdep Epolar & + &Epair Etriplet Equartet' + + + to_ev_per_atom = lo_Hartree_to_eV / ss%na + + do i = 1, sim%nt + total_energy = sum(pebuf(i,:)) + U0 + write (u, "(1X,I5,6(2X,E20.12))") i, sim%stat%potential_energy(i)*to_ev_per_atom, total_energy*to_ev_per_atom, pebuf(i, 4)*to_ev_per_atom, & + pebuf(i, 1)*to_ev_per_atom, pebuf(i, 2)*to_ev_per_atom, & + pebuf(i, 3)*to_ev_per_atom + end do + + ! close outfile.energies + close (u) + + write (*, '(A)') ' ... energies writen to `outfile.energies`' + end if -end block io +end block epotthings call mw%destroy() diff --git a/src/effective_hamiltonian/options.f90 b/src/effective_hamiltonian/options.f90 index 93a196dd..d8b52550 100644 --- a/src/effective_hamiltonian/options.f90 +++ b/src/effective_hamiltonian/options.f90 @@ -37,17 +37,22 @@ subroutine parse(opts) version=lo_version, & license=lo_licence, & help='Usage: ', & - description='Calculates the TDEP effective Hamiltonian from force constants. By default only the potential energy component is calculated.', & - examples=["mpirun effective_hamiltonian --thirdorder --fourthorder"], & + 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 anharmonic correction to the free energy', & + 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 anharmonic correction to the free energy', & + 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', & diff --git a/src/libolle/lo_epot.f90 b/src/libolle/lo_epot.f90 index ef5f154e..f7bc80c5 100644 --- a/src/libolle/lo_epot.f90 +++ b/src/libolle/lo_epot.f90 @@ -91,10 +91,10 @@ subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, m ! Calculate the energy call pot%energies_and_forces(p%u, e2, e3, e4, ep, f2, f3, f4, fp) - ebuf(i, 1) = e2 / ss%na - ebuf(i, 2) = e3 / ss%na - ebuf(i, 3) = e4 / ss%na - ebuf(i, 4) = ep / ss%na + ebuf(i, 1) = e2 + ebuf(i, 2) = e3 + ebuf(i, 3) = e4 + ebuf(i, 4) = ep end do call mem%deallocate(f2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) From 728bee6290e49a976ac4717d673b0353d4bb322c Mon Sep 17 00:00:00 2001 From: ejmeitz Date: Fri, 28 Mar 2025 20:41:26 -0400 Subject: [PATCH 4/8] got lost in my massive git mistake --- src/anharmonic_free_energy/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/anharmonic_free_energy/Makefile b/src/anharmonic_free_energy/Makefile index 2d288a40..079460ad 100644 --- a/src/anharmonic_free_energy/Makefile +++ b/src/anharmonic_free_energy/Makefile @@ -24,7 +24,7 @@ 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 $@ From bcfab2df1f0ba73d49ae6e2df07d095ad4ec91a4 Mon Sep 17 00:00:00 2001 From: ejmeitz Date: Fri, 28 Mar 2025 20:54:10 -0400 Subject: [PATCH 5/8] fix CI build --- src/effective_hamiltonian/Makefile | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/effective_hamiltonian/Makefile b/src/effective_hamiltonian/Makefile index f085e247..9eb9b4e8 100644 --- a/src/effective_hamiltonian/Makefile +++ b/src/effective_hamiltonian/Makefile @@ -7,7 +7,7 @@ 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 = -lolle -lflap $(blaslapackLIBS) $(incLIBShdf) $(incLIBSmpi) +LIBS = ../../lib/libolle.a -lflap $(blaslapackLIBS) $(incLIBShdf) $(incLIBSmpi) #OPT = -O0 -fbacktrace -fcheck=all -Wall F90 = $(FC) $(LPATH) $(IPATH) $(MODULE_FLAG) $(OBJECT_PATH) @@ -25,7 +25,7 @@ clean: $(F90) $(F90FLAGS) -c $< $(LIBS) $(OBJECT_PATH)main.o: main.f90 $(OBJECT_PATH)options.o - $(F90) $(F90FLAGS) -c main.f90 $(LIBS) -o $@ + $(F90) $(F90FLAGS) -c main.f90 -o $@ $(OBJECT_PATH)options.o: - $(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@ + $(F90) $(F90FLAGS) -c options.f90 -o $@ From 3f5238c6a5bf9ebdcbed91a910ab798b48f14c11 Mon Sep 17 00:00:00 2001 From: ejmeitz Date: Mon, 31 Mar 2025 11:44:03 -0400 Subject: [PATCH 6/8] fix harmonic average --- src/effective_hamiltonian/main.f90 | 213 +++++++++++++++++------------ src/libolle/lo_epot.f90 | 23 ++-- 2 files changed, 138 insertions(+), 98 deletions(-) diff --git a/src/effective_hamiltonian/main.f90 b/src/effective_hamiltonian/main.f90 index 7cb9e3a4..63095c2b 100644 --- a/src/effective_hamiltonian/main.f90 +++ b/src/effective_hamiltonian/main.f90 @@ -36,14 +36,19 @@ program effective_hamiltonian 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 - ! Get CLI options - call opts%parse() + + if (.not. mw%talk) opts%verbosity = -100 if (opts%nconf .gt. 0) generate_configs = .true. @@ -53,15 +58,13 @@ program effective_hamiltonian write (*, '(1X,A40,L3)') 'Fourthorder contribution ', opts%fourthorder write (*, '(1X,A40,I5)') 'Stride ', opts%stride write(*, '(1X,A40,L3)') 'Generate canonical configs ', generate_configs - write (*, '(1X,A40,L3)') 'Quantum configurations ', opts%quantum - write (*, '(1X,A40,F20.12)') 'Temperature ', opts%temperature + if(generate_configs) then + write(*, '(1X,A40,I5)') 'Num Configs ', opts%nconf + write (*, '(1X,A40,L3)') 'Quantum configurations ', opts%quantum + write (*, '(1X,A40,F20.12)') 'Temperature ', opts%temperature + end if end if - - call mw%init() - call mem%init() - if (.not. mw%talk) opts%verbosity = -100 - ! Read structures if (mw%talk) write (*, *) '... reading infiles' call ss%readfromfile('infile.ssposcar', verbosity=opts%verbosity) @@ -82,6 +85,11 @@ program effective_hamiltonian 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) @@ -105,24 +113,25 @@ program effective_hamiltonian energy : block - integer :: i, u + integer :: i,j, 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__) - ! call mem%allocate(kebuf, [opts%nconf], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) pebuf = 0.0_r8 - ! kebuf = 0.0_r8 - call pot%statistical_sampling(uc, ss, fc2, opts%nconf, opts%temperature, pebuf, mw, mem, opts%verbosity) + 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__) - ! call mem%allocate(kebuf, [sim%nt], 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 @@ -148,9 +157,10 @@ program effective_hamiltonian pebuf(i, 3) = e4 pebuf(i, 4) = ep end do - end if - call mw%allreduce('sum', pebuf) + call mw%allreduce('sum', pebuf) + + end if if (mw%talk) write (*, *) '... calculated energies' @@ -169,89 +179,116 @@ program effective_hamiltonian T_actual = sim%temperature_thermostat end if - ! 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 + ! Don't have sim%stat%potential energy when generating canonical_configs + ! So we cannot get U0 or cumulants + if (.not. generate_configs) then + ! Calculate the baseline energy + allocate (ediff(sim%nt, 5)) + ediff = 0.0_r8 - ! 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/real(ss%na, r8) - 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)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 1)*lo_Hartree_to_eV - write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2): ', cumulant(1, 2)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 2)*lo_Hartree_to_eV - write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar): ', cumulant(1, 3)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 3)*lo_Hartree_to_eV - if (opts%thirdorder) then - write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar-E3): ', cumulant(1, 4)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 4)*lo_Hartree_to_eV + 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 - if (opts%fourthorder) then - write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar-E3-E4): ', cumulant(1, 5)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 5)*lo_Hartree_to_eV + + ! 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/real(ss%na, r8) + 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)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 1)*lo_Hartree_to_eV + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2): ', cumulant(1, 2)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 2)*lo_Hartree_to_eV + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar): ', cumulant(1, 3)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 3)*lo_Hartree_to_eV + if (opts%thirdorder) then + write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar-E3): ', cumulant(1, 4)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 4)*lo_Hartree_to_eV + 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)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 5)*lo_Hartree_to_eV + end if + + close(u) + write (*, '(A)') ' ... cumulants writen to `outfile.cumulants`' 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)*lo_Hartree_to_eV + else if(opts%fourthorder) then + U0 = cumulant(1,5)*lo_Hartree_to_eV + else + U0 = cumulant(1,3)*lo_Hartree_to_eV + end if - if(opts%thirdorder .and. (.not. opts%fourthorder)) then - U0 = cumulant(1,4)*lo_Hartree_to_eV - else if(opts%fourthorder) then - U0 = cumulant(1,5)*lo_Hartree_to_eV - else - U0 = cumulant(1,3)*lo_Hartree_to_eV - 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: ', 'eV/atom' + write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) + write (u, *) 'U0 [eV/atom]: ', U0 + write (u, "(A)") '# conf Etotal_actual Etotal_tdep Epolar & + &Epair Etriplet Equartet' + + + to_ev_per_atom = lo_Hartree_to_eV / ss%na + + do i = 1, sim%nt + total_energy = sum(pebuf(i,:)) + U0 + write (u, "(1X,I5,6(2X,E20.12))") i, sim%stat%potential_energy(i)*to_ev_per_atom, total_energy*to_ev_per_atom, pebuf(i, 4)*to_ev_per_atom, & + pebuf(i, 1)*to_ev_per_atom, pebuf(i, 2)*to_ev_per_atom, & + pebuf(i, 3)*to_ev_per_atom + end do + + ! close outfile.energies + close (u) + + write (*, '(A)') ' ... energies writen to `outfile.energies`' + 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: ', 'eV/atom' - write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) - write (u, *) 'U0 [eV/atom]: ', U0 - write (u, "(A)") '# conf Etotal_actual Etotal_tdep Epolar & - &Epair Etriplet Equartet' + else + if (mw%talk) then + u = open_file('out', 'outfile.energies') + write (u, '(A,A)') '# Unit: ', 'eV/atom' + 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' - to_ev_per_atom = lo_Hartree_to_eV / ss%na - - do i = 1, sim%nt - total_energy = sum(pebuf(i,:)) + U0 - write (u, "(1X,I5,6(2X,E20.12))") i, sim%stat%potential_energy(i)*to_ev_per_atom, total_energy*to_ev_per_atom, pebuf(i, 4)*to_ev_per_atom, & - pebuf(i, 1)*to_ev_per_atom, pebuf(i, 2)*to_ev_per_atom, & - pebuf(i, 3)*to_ev_per_atom - end do + to_ev_per_atom = lo_Hartree_to_eV / ss%na + + do i = 1, opts%nconf + write (u, "(1X,I5,4(2X,E20.12))") i, pebuf(i, 4)*to_ev_per_atom, pebuf(i, 1)*to_ev_per_atom, & + pebuf(i, 2)*to_ev_per_atom, pebuf(i, 3)*to_ev_per_atom + end do - ! close outfile.energies - close (u) + close (u) - write (*, '(A)') ' ... energies writen to `outfile.energies`' + write (*, '(A)') ' ... energies writen to `outfile.energies`' + end if end if + end block epotthings diff --git a/src/libolle/lo_epot.f90 b/src/libolle/lo_epot.f90 index f7bc80c5..698fc2e0 100644 --- a/src/libolle/lo_epot.f90 +++ b/src/libolle/lo_epot.f90 @@ -33,7 +33,7 @@ module lo_epot contains !> statistically sample - subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, mem, verbosity) + 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 @@ -46,6 +46,10 @@ subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, m 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 @@ -55,15 +59,12 @@ subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, m type(lo_crystalstructure) :: p integer :: ctr, i - real(r8), dimension(:, :), allocatable, intent(out) :: ebuf real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp - real(r8) :: e2, e3, e4, ep, tomev, emean + real(r8) :: e2, e3, e4, ep ! Copy of structure to work with p = ss - ! Some helper space - call mem%allocate(ebuf, [nstep, 4], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) ebuf = 0.0_r8 ! Dummy space for force @@ -76,26 +77,28 @@ subroutine statistical_sampling(pot, uc, ss, fc, nstep, temperature, ebuf, mw, m f4 = 0.0_r8 fp = 0.0_r8 - ctr = 0 do i = 1, nstep - ctr = ctr + 1 ! does parallelizing like this create a copy of the IFCs on each rank? not ideal - if (mod(ctr, mw%n) .ne. mw%r) cycle + 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, .true., .false., -1.0_r8, mw, nosync=.true.) + 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 + 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__) From 40352d57f57ff45bb536c17954f7650d4a6d012d Mon Sep 17 00:00:00 2001 From: ejmeitz Date: Mon, 31 Mar 2025 12:52:08 -0400 Subject: [PATCH 7/8] fix some units stuff --- src/effective_hamiltonian/main.f90 | 53 ++++++++++++++---------------- 1 file changed, 25 insertions(+), 28 deletions(-) diff --git a/src/effective_hamiltonian/main.f90 b/src/effective_hamiltonian/main.f90 index 63095c2b..4197c249 100644 --- a/src/effective_hamiltonian/main.f90 +++ b/src/effective_hamiltonian/main.f90 @@ -113,7 +113,7 @@ program effective_hamiltonian energy : block - integer :: i,j, u + integer :: i, u real(r8), dimension(:, :), allocatable :: f2, f3, f4, fp real(r8) :: e2, e3, e4, ep, total_energy, to_ev_per_atom @@ -170,9 +170,11 @@ program effective_hamiltonian epotthings: block real(r8), dimension(3, 5) :: cumulant real(r8), dimension(:, :), allocatable :: ediff - real(r8) :: inverse_kbt, T_actual, U0, total_energy, to_ev_per_atom + 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 @@ -212,21 +214,21 @@ program effective_hamiltonian end do ! And normalize it to be per atom - cumulant = cumulant/real(ss%na, r8) + cumulant = cumulant*to_mev_per_atom if (mw%talk) then u = open_file('out', 'outfile.cumulants') - write (u, *) 'Temperature (K) : ', T_actual + 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)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 1)*lo_Hartree_to_eV - write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2): ', cumulant(1, 2)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 2)*lo_Hartree_to_eV - write (u, "(1X,A,E21.14,1X,A,F21.14)") ' mean(E-E2-Epolar): ', cumulant(1, 3)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 3)*lo_Hartree_to_eV + 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)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 4)*lo_Hartree_to_eV + 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)*lo_Hartree_to_eV, ' upper bound:', cumulant(2, 5)*lo_Hartree_to_eV + 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) @@ -234,30 +236,27 @@ program effective_hamiltonian end if if(opts%thirdorder .and. (.not. opts%fourthorder)) then - U0 = cumulant(1,4)*lo_Hartree_to_eV + U0 = cumulant(1,4) else if(opts%fourthorder) then - U0 = cumulant(1,5)*lo_Hartree_to_eV + U0 = cumulant(1,5) else - U0 = cumulant(1,3)*lo_Hartree_to_eV + 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: ', 'eV/atom' + write (u, '(A,A)') '# Unit: ', 'meV/atom' write (u, '(A,A)') '# no. atoms: ', tochar(ss%na) - write (u, *) 'U0 [eV/atom]: ', U0 + write (u, *) '# U0 [meV/atom]: ', U0 write (u, "(A)") '# conf Etotal_actual Etotal_tdep Epolar & &Epair Etriplet Equartet' - - - to_ev_per_atom = lo_Hartree_to_eV / ss%na do i = 1, sim%nt - total_energy = sum(pebuf(i,:)) + U0 - write (u, "(1X,I5,6(2X,E20.12))") i, sim%stat%potential_energy(i)*to_ev_per_atom, total_energy*to_ev_per_atom, pebuf(i, 4)*to_ev_per_atom, & - pebuf(i, 1)*to_ev_per_atom, pebuf(i, 2)*to_ev_per_atom, & - pebuf(i, 3)*to_ev_per_atom + total_energy = sum(pebuf(i,:))*to_mev_per_atom + U0 + write (u, "(1X,I5,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 @@ -269,18 +268,16 @@ program effective_hamiltonian else if (mw%talk) then u = open_file('out', 'outfile.energies') - write (u, '(A,A)') '# Unit: ', 'eV/atom' + 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' - - - to_ev_per_atom = lo_Hartree_to_eV / ss%na do i = 1, opts%nconf - write (u, "(1X,I5,4(2X,E20.12))") i, pebuf(i, 4)*to_ev_per_atom, pebuf(i, 1)*to_ev_per_atom, & - pebuf(i, 2)*to_ev_per_atom, pebuf(i, 3)*to_ev_per_atom + write (u, "(1X,I5,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) From 53c390a858607552f3531c7ce63c4e08295c2296 Mon Sep 17 00:00:00 2001 From: ejmeitz Date: Thu, 3 Apr 2025 10:07:40 -0400 Subject: [PATCH 8/8] print I8 --- src/effective_hamiltonian/main.f90 | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/effective_hamiltonian/main.f90 b/src/effective_hamiltonian/main.f90 index 4197c249..7ff1f8a4 100644 --- a/src/effective_hamiltonian/main.f90 +++ b/src/effective_hamiltonian/main.f90 @@ -59,7 +59,7 @@ program effective_hamiltonian write (*, '(1X,A40,I5)') 'Stride ', opts%stride write(*, '(1X,A40,L3)') 'Generate canonical configs ', generate_configs if(generate_configs) then - write(*, '(1X,A40,I5)') 'Num Configs ', opts%nconf + 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 @@ -181,8 +181,7 @@ program effective_hamiltonian T_actual = sim%temperature_thermostat end if - ! Don't have sim%stat%potential energy when generating canonical_configs - ! So we cannot get U0 or cumulants + if (.not. generate_configs) then ! Calculate the baseline energy allocate (ediff(sim%nt, 5)) @@ -254,7 +253,7 @@ program effective_hamiltonian do i = 1, sim%nt total_energy = sum(pebuf(i,:))*to_mev_per_atom + U0 - write (u, "(1X,I5,6(2X,E20.12))") i, sim%stat%potential_energy(i)*to_mev_per_atom, total_energy, pebuf(i, 4)*to_mev_per_atom, & + 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 @@ -264,7 +263,9 @@ program effective_hamiltonian 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') @@ -276,7 +277,7 @@ program effective_hamiltonian &Epair Etriplet Equartet' do i = 1, opts%nconf - write (u, "(1X,I5,4(2X,E20.12))") i, pebuf(i, 4)*to_mev_per_atom, pebuf(i, 1)*to_mev_per_atom, & + 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