#include "precompilerdefinitions" program extract_forceconstants !!{!src/extract_forceconstants/manual.md!}! use konstanter, only: flyt,lo_tol,lo_bohr_to_A,lo_pressure_HartreeBohr_to_GPa,lo_Hartree_to_eV use helpers, only: lo_mpi_helper,tochar,walltime,open_file,lo_chop,lo_does_file_exist,lo_trueNtimes,lo_mean,lo_stddev use type_crystalstructure, only: lo_crystalstructure use type_forceconstant_firstorder, only: lo_forceconstant_firstorder 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_jij_secondorder, only: lo_jij_secondorder use type_symmetrylist, only: lo_symlist use type_forcemap, only: lo_forcemap,lo_secondorder_rot_herm_huang use type_mdsim, only: lo_mdsim ! use options, only: lo_opts use test_forceconstant_symmetry, only: test_forceconstant_constraints use relax_structure, only: relax_internal_coordinates use dipole, only: subtract_longrange_interactions,histfit ! implicit none ! type(lo_opts) :: opts type(lo_crystalstructure) :: uc,ss type(lo_mdsim) :: sim type(lo_forcemap) :: map type(lo_forceconstant_firstorder) :: fc1 type(lo_forceconstant_secondorder) :: fc2 type(lo_forceconstant_thirdorder) :: fc3 type(lo_forceconstant_fourthorder) :: fc4 type(lo_jij_secondorder) :: jij type(lo_mpi_helper) :: mw real(flyt) :: t0 ! I should determine the general nature of the task here: it could be ! a) ! get forceconstants from force-displacement data ! b) ! get forceconstants from an input file, for interpolation purposes ! To start, I need at least the unitcell call mw%init() t0=walltime() call opts%parse() write(*,*) '... reading unitcell' call uc%readfromfile('infile.ucposcar',verbosity=opts%verbosity) call uc%classify('wedge',timereversal=.true.) ! We need a forcemap to continue. getfrcmap: block type(lo_symlist) :: sl integer :: i real(flyt) :: t0 ! Either read it from file or create it if ( opts%readforcemap ) then ! I might need the supercell anyway if ( opts%readirreducible .eqv. .false. ) then call ss%readfromfile('infile.ssposcar') call ss%classify('supercell',uc) endif call map%read_from_hdf5(uc,ss,'infile.forcemap.hdf5',opts%readirreducible,opts%verbosity) ! Update constraints that depend on structure call lo_secondorder_rot_herm_huang(map,uc,map%constraints%eq2,map%constraints%neq2,& opts%rotationalconstraints,opts%huanginvariance,opts%hermitian) else ! Now we have to calculate the whole thing. call uc%classify('wedge',timereversal=.true.) call ss%readfromfile('infile.ssposcar') write(*,*) '... min cutoff: ',tochar(ss%mincutoff()*lo_bohr_to_A) write(*,*) '... max cutoff: ',tochar(ss%maxcutoff()*lo_bohr_to_A) ! Get all the symmetries call sl%generate(uc,ss,opts%cutoff2,opts%cutoff3,opts%cutoff4,opts%polar,& transposition=opts%transposition,spacegroup=opts%spacegroup,verbosity=opts%verbosity+1,& wzdim=opts%wzdimensions,nj2=opts%njump2,nj3=opts%njump3,& nj4=opts%njump4,firstorder=opts%firstorder,magcutoff2=opts%magcutoff2,magsinglet=opts%magnetic_onsite) ! And create the map call map%generate(sl,uc,ss,polarcorrectiontype=opts%polarcorrectiontype,verbosity=opts%verbosity) ! Create the constraints t0=walltime() write(*,*) '... generating constraints' t0=walltime() call map%forceconstant_constraints(uc,opts%rotationalconstraints,opts%huanginvariance,opts%hermitian,opts%verbosity+10) write(*,*) '... got constraints in ',tochar(walltime()-t0),'s' if ( opts%printforcemap ) then call map%write_to_hdf5(sl,uc,ss,'outfile.forcemap.hdf5',opts%verbosity) endif endif ! Now that we have the map, create the linear constraints end block getfrcmap ! With the forcemap, we can get the irreducible representation getirr: block real(flyt) :: sz integer :: i,u ! Either read them from file, or calculate them if ( opts%readirreducible ) then if ( map%firstorder ) then u=open_file('in','infile.irrifc_firstorder') do i=1,map%nfc_singlet read(u,*) map%ifc_singlet(i) enddo close(u) endif if ( map%secondorder ) then u=open_file('in','infile.irrifc_secondorder') do i=1,map%nfc_pair read(u,*) map%ifc_pair(i) enddo close(u) endif if ( map%thirdorder ) then u=open_file('in','infile.irrifc_thirdorder') do i=1,map%nfc_triplet read(u,*) map%ifc_triplet(i) enddo close(u) endif if ( map%fourthorder ) then u=open_file('in','infile.irrifc_fourthorder') do i=1,map%nfc_quartet read(u,*) map%ifc_quartet(i) enddo close(u) endif if ( map%polar ) then u=open_file('in','infile.irrtheta_loto') do i=1,map%ntheta_eps read(u,*) map%theta_eps(i) enddo do i=1,map%ntheta_Z read(u,*) map%theta_Z(i) enddo close(u) endif if ( map%magnetic_pair_interactions ) then write(*,*) 'FIXME READ IRREDUCIBLE MAGNETIC PAIRS' stop endif write(*,*) '... enforce the constraints' call map%enforce_constraints() else ! Now we have to calculate the whole thing, first we need positions and forces. ! Get the force-displacement data if ( lo_does_file_exist('infile.sim.hdf5') ) then call sim%read_from_hdf5('infile.sim.hdf5',verbosity=2,stride=opts%stride) else if ( map%magnetic_pair_interactions ) then call sim%read_from_file(verbosity=2,stride=opts%stride,readmag=.true.) else call sim%read_from_file(verbosity=2,stride=opts%stride,readmag=.false.) endif call sim%rms_displacement(uc,ss) endif ! Fix drift call sim%remove_force_and_center_of_mass_drift() !call histfit(sim) ! Solve for the born-charges and dielectric tensor first. if ( map%polar ) then ! get the raw Born charges and dielectric tensor call map%solve_for_borncharges(filename='infile.lotosplitting',verbosity=opts%verbosity+1) endif ! Remove the dipole-dipole interactions if ( map%polarcorrectiontype .eq. 3 ) then call subtract_longrange_interactions(sim,map,uc,ss) endif ! Now I can figure out how much memory is needed for this, approximately. It's nice to ! warn the use here if it goes into TB range or something. sz=(sim%na*3*sim%nt*map%nfc*4*3)/1024_flyt**2 if ( sz .gt. 1000 ) then write(*,*) '... will try to allocate at least ',tochar(ceiling(sz)),'Mb' endif ! Maybe relax? if ( opts%relax ) then call relax_internal_coordinates(map,sim,uc,ss,opts%solver,opts%solveralpha) ! Dump the relaxed structures call uc%writetofile('outfile.relaxed_ucposcar',1) call ss%writetofile('outfile.relaxed_ssposcar',1) write(*,*) 'Relaxation done, files written!' stop endif ! Solve for the irreducible representation call map%solve_equations(sim,opts%solver,opts%solveralpha,opts%residualfit,verbosity=2) ! Solve for magnetic Jij if ( map%magnetic_pair_interactions ) then call map%solve_for_jij(sim,uc,ss,verbosity=2) endif ! And store the irreducible representation if ( map%firstorder ) then u=open_file('out','outfile.irrifc_firstorder') do i=1,map%nfc_singlet write(u,*) map%ifc_singlet(i) enddo close(u) endif if ( map%secondorder ) then u=open_file('out','outfile.irrifc_secondorder') do i=1,map%nfc_pair write(u,*) map%ifc_pair(i) enddo close(u) endif if ( map%thirdorder ) then u=open_file('out','outfile.irrifc_thirdorder') do i=1,map%nfc_triplet write(u,*) map%ifc_triplet(i) enddo close(u) endif if ( map%fourthorder ) then u=open_file('out','outfile.irrifc_fourthorder') do i=1,map%nfc_quartet write(u,*) map%ifc_quartet(i) enddo close(u) endif if ( map%polar ) then u=open_file('out','outfile.irrtheta_loto') do i=1,map%ntheta_eps write(u,*) map%theta_eps(i) enddo do i=1,map%ntheta_Z write(u,*) map%theta_Z(i) enddo close(u) endif if ( map%magnetic_pair_interactions ) then u=open_file('out','outfile.irrtheta_magsecondorder') do i=1,map%ntheta_magpair_jij write(u,*) map%theta_magpair_jij(i) enddo do i=1,map%ntheta_magpair_tij write(u,*) map%theta_magpair_tij(i) enddo close(u) endif endif end block getirr ! Now we can get the actual forceconstants from the forcemap and the irreducible representation getfc: block integer :: i if ( map%firstorder ) then call map%get_firstorder_forceconstant(uc,fc1) write(*,*) 'residual forces (eV/A):' do i=1,uc%na write(*,*) i,fc1%atom(i)%m enddo endif if ( map%secondorder ) then call map%get_secondorder_forceconstant(uc,fc2) call fc2%writetofile(uc,'outfile.forceconstant') call fc2%get_elastic_constants(uc) write(*,*) 'elastic constants (GPa):' do i=1,6 write(*,"(6(3X,F15.5))") lo_chop(fc2%elastic_constants_voigt(:,i)*lo_pressure_HartreeBohr_to_GPa,lo_tol) enddo endif if ( map%thirdorder ) then call map%get_thirdorder_forceconstant(uc,fc3) call fc3%writetofile(uc,'outfile.forceconstant_thirdorder') endif if ( map%fourthorder ) then call map%get_fourthorder_forceconstant(uc,fc4) call fc4%writetofile(uc,'outfile.forceconstant_fourthorder') endif if ( map%magnetic_pair_interactions ) then call map%get_secondorder_jij(uc,jij) call jij%writetofile(uc,'outfile.jij') endif end block getfc ! And perhaps calculate U0 if ( opts%ediff ) then getU0: block type(lo_forceconstant_secondorder) :: fc2_ss type(lo_forceconstant_thirdorder) :: fc3_ss type(lo_forceconstant_fourthorder) :: fc4_ss real(flyt), dimension(:,:,:,:), allocatable :: polarfc real(flyt), dimension(:,:), allocatable :: fp,f2,f3,f4,de,dle 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 :: t,i,a1,a2,a3,a4,i1,i2,i3,i4,u ! Make sure we have the supercell and MD data if ( ss%na .lt. 0 ) then call ss%readfromfile('infile.ssposcar') call ss%classify('supercell',uc) endif if ( sim%na .lt. 0 ) then call sim%read_from_file(verbosity=2,stride=opts%stride) call sim%remove_force_and_center_of_mass_drift() endif ! Get forceconstants for the supercell if ( map%secondorder ) call fc2%remap(uc,ss,fc2_ss) if ( map%thirdorder ) call fc3%remap(uc,ss,fc3_ss) if ( map%fourthorder ) call fc4%remap(uc,ss,fc4_ss) if ( map%polar ) then lo_allocate(polarfc(3,3,ss%na,ss%na)) call fc2%supercell_longrange_dynamical_matrix_at_gamma(ss,polarfc,1E-15_flyt) endif ! Space for forces and stuff lo_allocate(fp(3,ss%na)) lo_allocate(f2(3,ss%na)) lo_allocate(f3(3,ss%na)) lo_allocate(f4(3,ss%na)) lo_allocate(de(sim%nt,5)) lo_allocate(dle(sim%nt,5)) fp=0.0_flyt f2=0.0_flyt f3=0.0_flyt f4=0.0_flyt de=0.0_flyt dle=0.0_flyt write(*,*) '' write(*,*) 'CALCULATING POTENTIAL ENERGIES (eV/atom)' write(*,*) ' conf Epot Epolar Epair Etriplet Equartet' ! 1 -6.839524677083 0.016417795628 0.112660623299 0.001685955875 0.002296712482 ! 2 -6.829915786458 0.019703035642 0.118811923784 0.003986464042 0.001081803009 ! Calculate energies and stuff do t=1,sim%nt ! get the raw energy de(t,1)=sim%stat%ep(t) ! secondorder energy, first the polar term (if applicable) if ( fc2%polar ) then energy=0.0_flyt do a1=1,ss%na v0=0.0_flyt do a2=1,ss%na v0=v0+matmul(polarfc(:,:,a2,a1),sim%u(:,a2,t)) enddo energy=energy+dot_product(sim%u(:,a1,t),v0)*0.5_flyt enddo de(t,2)=energy endif ! then the pair term if ( map%secondorder ) then energy=0.0_flyt do a1=1,ss%na v0=0.0_flyt do i1=1,fc2_ss%atom(a1)%n a2=fc2_ss%atom(a1)%pair(i1)%i2 m2=fc2_ss%atom(a1)%pair(i1)%m v0=v0+matmul(m2,sim%u(:,a2,t)) enddo energy=energy+dot_product(sim%u(:,a1,t),v0)*0.5_flyt enddo de(t,3)=energy endif ! triplet term if ( map%thirdorder ) then f3=0.0_flyt energy=0.0_flyt do a1=1,fc3_ss%na v0=0.0_flyt do i=1,fc3_ss%atom(a1)%n m3=fc3_ss%atom(a1)%triplet(i)%m a2=fc3_ss%atom(a1)%triplet(i)%i2 a3=fc3_ss%atom(a1)%triplet(i)%i3 u2=sim%u(:,a2,t) u3=sim%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),sim%u(:,a1,t))/3.0_flyt enddo de(t,4)=energy endif ! quartet term if ( map%fourthorder ) then f4=0.0_flyt energy=0.0_flyt do a1=1,fc4_ss%na v0=0.0_flyt do i=1,fc4_ss%atom(a1)%n m4=fc4_ss%atom(a1)%quartet(i)%m a2=fc4_ss%atom(a1)%quartet(i)%i2 a3=fc4_ss%atom(a1)%quartet(i)%i3 a4=fc4_ss%atom(a1)%quartet(i)%i4 u2=sim%u(:,a2,t) u3=sim%u(:,a3,t) u4=sim%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),sim%u(:,a1,t))/4.0_flyt enddo de(t,5)=energy endif if ( lo_trueNtimes(t,20,sim%nt) ) then write(*,"(1X,I5,5(2X,F20.12))") t,de(t,:)*lo_Hartree_to_eV/ss%na endif enddo ! Calculate some differences dle(:,1)=de(:,1) dle(:,2)=de(:,1)-de(:,2) dle(:,3)=de(:,1)-de(:,2)-de(:,3) dle(:,4)=de(:,1)-de(:,2)-de(:,3)-de(:,4) dle(:,5)=de(:,1)-de(:,2)-de(:,3)-de(:,4)-de(:,5) ! Convert to eV/atom dle=dle*lo_Hartree_to_eV/sim%na write(*,*) '' write(*,"(1X,' avg E',F20.12,' (eV/atom) stddev ',F20.12,' eV/atom')") lo_mean(dle(:,1)),lo_stddev(dle(:,1)) write(*,"(1X,' avg E-Epolar',F20.12,' (eV/atom) stddev ',F20.12,' eV/atom')") lo_mean(dle(:,2)),lo_stddev(dle(:,2)) write(*,"(1X,' avg E-Epolar-Epair',F20.12,' (eV/atom) stddev ',F20.12,' eV/atom')") lo_mean(dle(:,3)),lo_stddev(dle(:,3)) write(*,"(1X,' avg E-Epolar-Epair-Etriplet',F20.12,' (eV/atom) stddev ',F20.12,' eV/atom')") lo_mean(dle(:,4)),lo_stddev(dle(:,4)) write(*,"(1X,' avg E-Epolar-Epair-Etriplet-Equartet',F20.12,' (eV/atom) stddev ',F20.12,' eV/atom')") lo_mean(dle(:,5)),lo_stddev(dle(:,5)) ! Dump it to file u=open_file('out','outfile.U0') write(u,"(4(1X,E19.12))") lo_mean(dle(:,1)),lo_mean(dle(:,3)),lo_mean(dle(:,4)),lo_mean(dle(:,5)) close(u) end block getU0 endif ! Now, perhaps test all the symmetry stuff on the forceconstants !call test_forceconstant_constraints(map,fc1,fc2,fc3,fc4,uc) write(*,*) 'Done in ',tochar(walltime()-t0),'s' call mw%destroy() end program