From 29f3312cfc6c3de05c4788bc8015cfb0c7e077df Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 31 Jul 2023 14:47:09 +0200 Subject: [PATCH] remove useless routines --- input/methods | 2 +- input/options | 2 +- src/GT/G0T0pp.f90 | 15 +- src/GT/GTpp_self_energy_diag.f90 | 1 + utils/install_lapack.sh | 7 - utils/reblock.f90 | 2764 ------------------------------ 6 files changed, 17 insertions(+), 2774 deletions(-) delete mode 100755 utils/install_lapack.sh delete mode 100644 utils/reblock.f90 diff --git a/input/methods b/input/methods index 1baccb4..91124e8 100644 --- a/input/methods +++ b/input/methods @@ -15,5 +15,5 @@ # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW F F F F F F # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh - F F F F F F + T F F F F F # * unrestricted version available diff --git a/input/options b/input/options index 7089880..0b1a49b 100644 --- a/input/options +++ b/input/options @@ -11,7 +11,7 @@ # GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg 256 0.00001 T 5 T 0.0 F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg - 256 0.00001 T 5 F 0.0 T F + 256 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index 9c382f1..c33dcbb 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -72,6 +72,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp ! Dimensions of the pp-RPA linear reponse matrices +! nOOs = nO*(nO + 1)/2 +! nVVs = nV*(nV + 1)/2 + nOOs = nO*nO nVVs = nV*nV @@ -93,6 +96,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp !---------------------------------------------- ispin = 1 +! iblock = 1 iblock = 3 ! Compute linear response @@ -115,6 +119,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp !---------------------------------------------- ispin = 2 +! iblock = 2 iblock = 4 ! Compute linear response @@ -133,15 +138,21 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp call print_excitation_energies('ppRPA (N-2)',iblock,nOOt,Om2t) !---------------------------------------------- -! Compute T-matrix version of the self-energy +! Compute excitation densities !---------------------------------------------- +! iblock = 1 iblock = 3 call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) +! iblock = 2 iblock = 4 call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) +!---------------------------------------------- +! Compute T-matrix version of the self-energy +!---------------------------------------------- + call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) @@ -173,6 +184,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp ! Compute the ppRPA correlation energy ispin = 1 +! iblock = 1 iblock = 3 allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) @@ -186,6 +198,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp deallocate(Bpp,Cpp,Dpp) ispin = 2 +! iblock = 2 iblock = 4 allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) diff --git a/src/GT/GTpp_self_energy_diag.f90 b/src/GT/GTpp_self_energy_diag.f90 index 3612d0e..9cd2a63 100644 --- a/src/GT/GTpp_self_energy_diag.f90 +++ b/src/GT/GTpp_self_energy_diag.f90 @@ -59,6 +59,7 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 enddo + enddo enddo diff --git a/utils/install_lapack.sh b/utils/install_lapack.sh deleted file mode 100755 index 203d7fb..0000000 --- a/utils/install_lapack.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/bash -x - -git clone https://github.com/Reference-LAPACK/lapack-release.git -cd lapack-release -cp make.inc.example make.inc -make -j 8 -mv librefblas.a liblapack.a libtmglib.a /Users/loos/Dropbox/quack/lib diff --git a/utils/reblock.f90 b/utils/reblock.f90 deleted file mode 100644 index 063344b..0000000 --- a/utils/reblock.f90 +++ /dev/null @@ -1,2764 +0,0 @@ -!--------------------------------------------------------------------! -! REBLOCK ! -! ======= ! -! Neil Drummond, 8.2005 ! -! (Based on earlier MDT reblock utility for old-style .hist files) ! -! ! -! This utility performs a statistical analysis of the raw QMC data ! -! held in the .hist file. The reblocking procedure is ! -! necessary in order to obtain reliable statistical error bars for ! -! the mean values of serially correlated data. Please refer to the ! -! CASINO manual for further information about the procedure. ! -! ! -! Changes ! -! NDD 9.05 Data files renamed vmc.hist, etc. Check for old format.! -! Other new checks on files. ! -! NDD 9.05 Changed format of .hist files (again). ! -! AB 5.06 Add reblocking analysis for future walking estimators. ! -! AB 11.07 Add reblocking analysis of forces and introduce a line ! -! break in the .hist file to allow storing many items. ! -! NDD 05.08 Rearranged output, to put important stuff at end. ! -! NDD 05.10 Allow for FISQ data in qmc.hist. ! -!--------------------------------------------------------------------! - - -MODULE stats_calcs -!-------------------------------------------------------------------! -! A collection of subroutines for performing various statistical ! -! analyses of data. ! -!-------------------------------------------------------------------! - IMPLICIT NONE - - -CONTAINS - - - SUBROUTINE compute_stats_unweighted(want_skew_kurt,n,data_arr,av,var,skew, & - &kurt,max_val,min_val) -!-------------------------------------------------------------------! -! Compute mean, variance, skewness, kurtosis and max and min of a ! -! set of data. ! -!-------------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: n - DOUBLE PRECISION,INTENT(in) :: data_arr(n) - LOGICAL,INTENT(in) :: want_skew_kurt - DOUBLE PRECISION,INTENT(out) :: av,var,skew,kurt,max_val,min_val - INTEGER i - DOUBLE PRECISION sum_delta_x2,sum_delta_x3,sum_delta_x4 - - if(n<2)then - write(6,*)'Can''t compute variance with fewer than two points.' - stop - endif - -! Compute average. - av=sum(data_arr(1:n))/dble(n) - -! Compute max and min. - max_val=maxval(data_arr(1:n)) - min_val=minval(data_arr(1:n)) - - if(want_skew_kurt)then -! Compute variance, skewness and kurtosis. - sum_delta_x2=0.d0 ; sum_delta_x3=0.d0 ; sum_delta_x4=0.d0 - do i=1,n - sum_delta_x2=sum_delta_x2+(data_arr(i)-av)**2 - sum_delta_x3=sum_delta_x3+(data_arr(i)-av)**3 - sum_delta_x4=sum_delta_x4+(data_arr(i)-av)**4 - enddo ! i - var=sum_delta_x2/dble(n-1) - if(var>0.d0)then - skew=((sqrt(dble(n-1))*dble(n))/dble(n-2))*sum_delta_x3/sum_delta_x2**1.5d0 - kurt=((dble(n+1)*dble(n)*dble(n-1))/(dble(n-2)*dble(n-3)))*& - &sum_delta_x4/sum_delta_x2**2-& - &((dble(n-1)*dble(n-1))/(dble(n-2)*dble(n-3)))*3.d0 - else - skew=0.d0 - kurt=0.d0 - endif - else -! Compute variance. - sum_delta_x2=0.d0 - do i=1,n - sum_delta_x2=sum_delta_x2+(data_arr(i)-av)**2 - enddo ! i - var=sum_delta_x2/dble(n-1) - skew=0.d0 - kurt=0.d0 - endif ! want_skew_kurt - - END SUBROUTINE compute_stats_unweighted - - - SUBROUTINE reblock_unweighted(no_pts,data_array,block_length,av, & - &std_err,delta_std_err) -!---------------------------------------------------------------! -! Compute the unweighted average of the data, and calculate the ! -! error bar for a given block length. ! -!---------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: no_pts,block_length - DOUBLE PRECISION,INTENT(in) :: data_array(:) - DOUBLE PRECISION,INTENT(out) :: av,std_err,delta_std_err - INTEGER i,k,no_blocks,no_pts_in_last_block,j - DOUBLE PRECISION last_block_weight,var,tot_weight, & - &tot_weight_sq,block_av,red_tot_weight,rec_block_length - -! Compute average of data. - av=sum(data_array(1:no_pts))/dble(no_pts) - -! Number of blocks. - no_blocks=no_pts/block_length - rec_block_length=1.d0/dble(block_length) - -! Evaluate the sum of the squares of the deviations from the average. -! Weight the last, incomplete block by its size as a fraction of the others. - var=0.d0 - k=0 - do i=1,no_blocks - block_av=0.d0 - do j=1,block_length - k=k+1 - block_av=block_av+data_array(k) - enddo ! j - block_av=block_av*rec_block_length - var=var+(block_av-av)**2 - enddo ! i - block_av=0.d0 - no_pts_in_last_block=0 - do - k=k+1 - if(k>no_pts)exit - no_pts_in_last_block=no_pts_in_last_block+1 - block_av=block_av+data_array(k) - enddo ! k - last_block_weight=dble(no_pts_in_last_block)*rec_block_length - if(no_pts_in_last_block>0)then - block_av=block_av/dble(no_pts_in_last_block) - var=var+(block_av-av)**2*last_block_weight - endif ! last block nonzero - -! Evaluate variance, standard error in mean and error in standard error. - tot_weight=dble(no_blocks)+last_block_weight - tot_weight_sq=dble(no_blocks)+last_block_weight**2 - red_tot_weight=tot_weight-tot_weight_sq/tot_weight - var=var/red_tot_weight - std_err=sqrt(var/tot_weight) - if(tot_weight>1.d0)then - delta_std_err=std_err/sqrt(2.d0*(tot_weight-1.d0)) - else - delta_std_err=0.d0 - endif - - END SUBROUTINE reblock_unweighted - - - SUBROUTINE reblock_weighted(no_pts,data_array,weight_array,block_length, & - &av,std_err,delta_std_err) -!--------------------------------------------------------------! -! Compute the weighted average of the data, and calculate the ! -! error bar for a given block length. ! -!--------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: no_pts,block_length - DOUBLE PRECISION,INTENT(in) :: weight_array(no_pts),data_array(no_pts) - DOUBLE PRECISION,INTENT(out) :: av,std_err,delta_std_err - INTEGER i,k,no_blocks,no_pts_in_last_block,j - DOUBLE PRECISION var,tot_weight, & - &tot_weight_sq,block_av,red_tot_weight,block_weight, & - &eff_no_blocks - -! Compute average of data. - av=0.d0 - tot_weight=0.d0 - do i=1,no_pts - av=av+data_array(i)*weight_array(i) - tot_weight=tot_weight+weight_array(i) - enddo ! i - av=av/tot_weight - -! Number of blocks - no_blocks=no_pts/block_length - -! Evaluate the sum of the squares of the deviations from the average. -! Last, incomplete block has fewer data points and hence a smaller weight. - var=0.d0 - tot_weight_sq=0.d0 - k=0 - do i=1,no_blocks - block_av=0.d0 - block_weight=0.d0 - do j=1,block_length - k=k+1 - block_av=block_av+data_array(k)*weight_array(k) - block_weight=block_weight+weight_array(k) - enddo ! j - block_av=block_av/block_weight - var=var+(block_av-av)**2*block_weight - tot_weight_sq=tot_weight_sq+block_weight**2 - enddo ! i - block_av=0.d0 - block_weight=0.d0 - no_pts_in_last_block=0 - do - k=k+1 - if(k>no_pts)exit - no_pts_in_last_block=no_pts_in_last_block+1 - block_av=block_av+data_array(k)*weight_array(k) - block_weight=block_weight+weight_array(k) - enddo ! k - if(no_pts_in_last_block>0)then - block_av=block_av/block_weight - var=var+(block_av-av)**2*block_weight - tot_weight_sq=tot_weight_sq+block_weight**2 - endif ! last block nonzero - -! Evaluate variance, standard error in mean and error in standard error. - red_tot_weight=tot_weight-tot_weight_sq/tot_weight - var=var/red_tot_weight - - eff_no_blocks=dble(no_blocks)+dble(no_pts_in_last_block)/dble(block_length) - - std_err=sqrt(var/eff_no_blocks) - if(eff_no_blocks>1.d0)then - delta_std_err=std_err/sqrt(2.d0*(eff_no_blocks-1.d0)) - else - delta_std_err=0.d0 - endif - - END SUBROUTINE reblock_weighted - - - SUBROUTINE correlation_time(n,Odata,Otau,Otau_err,Oave_in,Ovar_in) -!------------------------------------------------------------------------! -! Obtain correlation time from a set of data ! -!------------------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: n - DOUBLE PRECISION,INTENT(in) :: Odata(n) - DOUBLE PRECISION,INTENT(in),OPTIONAL :: Oave_in,Ovar_in - DOUBLE PRECISION,INTENT(out) :: Otau,Otau_err - DOUBLE PRECISION Oave,Oave2,O2ave,Ovar,invOvar,Oacorr,ri,invn - DOUBLE PRECISION,PARAMETER :: tol=1.d-100 - INTEGER i,sqrtn - - Otau=-1.d0 ; Otau_err=-1.d0 - if(n<10)return - invn=1.d0/dble(n) - sqrtn=nint(sqrt(dble(n))) - -! , **2, , variance - if(present(Oave_in))then - Oave=Oave_in - else - Oave=sum(Odata)*invn - endif - if(present(Ovar_in))then - Ovar=Ovar_in*invn*(n-1) - else - Oave2=Oave**2 - O2ave=sum(Odata**2)*invn - Ovar=O2ave-Oave2 - endif - if(Ovar tau - Otau=1.d0 - do i=1,n-1 - Oacorr=sum((Odata(1:n-i)-Oave)*(Odata(1+i:n)-Oave))*invOvar/dble(n-i) - Otau=Otau+2*Oacorr - if(i>=nint(3*Otau))then - ri=dble(i) ; exit - endif - enddo - -! Error in tau - Otau_err=Otau*sqrt((4*ri+2.d0)*invn) - - END SUBROUTINE correlation_time - - -END MODULE stats_calcs - - -MODULE analysis -!-------------------------------------------------------------! -! Miscellaneous subroutines for reading & analysing the data. ! -!-------------------------------------------------------------! - USE stats_calcs - IMPLICIT NONE - -! Tags for the columns of the data file, specifying where each data item -! is held. If a tag is negative, the data item isn't present. - INTEGER tag_step,tag_energy,tag_etotalt,tag_esqr,tag_popavgsqr,tag_K,tag_T, & - &tag_fisq,tag_Ewald,tag_local,tag_nonlocal,tag_short,tag_long,tag_cppei, & - &tag_cppe,tag_cppee,tag_masspol,tag_massvel,tag_darwinen,tag_darwinee, & - &tag_retard,tag_weight,tag_nconf,tag_eref,tag_ebest,tag_acc,tag_teff, & - &tag_dipole1,tag_dipole2,tag_dipole3,tag_dipole_sq,tag_contact_den, & - &tag_future0,tag_future1,tag_future2,tag_future3,tag_future4,tag_future5, & - &tag_future6,tag_future7,tag_future8,tag_future9,tag_future10 - -! Number of columns of data in .hist file. - INTEGER no_cols_qmc - -! Title of .hist file - CHARACTER(72) title - -! File version number - INTEGER version - -! CASINO input keywords: interaction type and basis type. - CHARACTER(20) interaction,atom_basis_type - -! Do we have Ewald/Coulomb interaction? Do we have MPC interaction? - LOGICAL coul_mpc,coul_ewald - -! Ion-ion energy - DOUBLE PRECISION constant_energy - -! Total number of electrons; no. of atoms per prim cell; no. primitive cells - INTEGER netot,nbasis,npcells - -! Number of parts in simulation cell: -! =npcells for periodic systems -! =netot for electron gas -! =1 otherwise - INTEGER nparts_per_simcell - -! Is the system periodic? - LOGICAL isperiodic - -! QMC method used - CHARACTER(3) qmc_method - -! Number of lines of data. - INTEGER Nlines - -! Number of equilibration lines. - INTEGER Nequil - -! Name of .hist file. - CHARACTER(8) filename - -! Number of initial lines to discard. - INTEGER Nskip - -! Array with the hist data from the files. - DOUBLE PRECISION,ALLOCATABLE :: data_array(:,:) - -! Energy units - CHARACTER(15) e_units - -! Units conversion: a.u.->eV and a.u.->kcal. - DOUBLE PRECISION,PARAMETER :: htoev=27.2113962d0,htokcal=627.507541278d0 - -! Are forces to be calculated? - INTEGER iion,iaxis,item,nitot_forces,naxis_forces,nitem_forces,& - &nitot_max_forces - INTEGER,ALLOCATABLE :: tag_forces(:,:,:) - DOUBLE PRECISION,ALLOCATABLE :: forces_array(:,:) - LOGICAL forces - - ! Use weights when calculating average energy, etc. - LOGICAL,PARAMETER :: use_weights=.true. - - -CONTAINS - - - SUBROUTINE read_header(io,dmc) -!----------------------------------------------------------------------! -! Read in the data in the .hist file header and count the lines, etc. ! -!----------------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: io - LOGICAL,INTENT(in) :: dmc - INTEGER ierr,i,s,isper_flag,nbreak,ialloc - CHARACTER(1) temp - CHARACTER(72) datastring - CHARACTER(500) checkline - LOGICAL,PARAMETER :: verbose=.false. - -! Check we don't have an old-style .hist file. - rewind(io) - read(io,'(a)',iostat=ierr)checkline - call check_ierr(ierr) - checkline=adjustl(checkline) - if(index(checkline,'Block')>0)then - write(6,*)'You appear to be analyzing a CASINO version 1 vmc.hist file. & - &Please use the' - write(6,*)'UPDATE_HIST utility to update it to the new format.' - stop - endif ! Old-style vmc.hist - if(index(checkline,'#')==0)then - write(6,*)'Your data file does not seem to start with a header. & - &This may be because you' - write(6,*)'are using an old-format file. If this is the case then & - &please use UPDATE_HIST' - write(6,*)'to update your file.' - stop - endif ! No header. - -! Count the data lines. Ignore comments. - rewind(io) - Nlines=0 - Nequil=0 - forces=.false. - do - read(io,'(a)',iostat=ierr)datastring - if(ierr>0)then - write(6,*)'Error reading data file.' - stop - endif - if(ierr<0)exit - if(index(datastring,'#')==0)then - Nlines=Nlines+1 - else - if(trim(adjustl(datastring))=='#### START STATS')Nequil=Nlines - endif ! Line not a comment. - if(.not.forces)then - if(index(datastring,'FOR')>0)forces=.true. ! atomic forces present - endif - enddo ! lines - if(dmc.and.Nequil==0)Nequil=Nlines - rewind(io) - if(verbose)then - if(Nlines/=1)then - write(6,*)'There are '//trim(i2s(Nlines))//' lines of data in ' & - &//trim(filename)//'.' - else - write(6,*)'There is 1 line of data in '//trim(filename)//'.' - endif ! Singular / plural - if(Nequil>1)then - write(6,*)'Of these, '//trim(i2s(Nlines))//' lines are marked as & - &equilibration data.' - elseif(Nequil==1)then - write(6,*)'Of these, 1 line is marked as equilibration data.' - else - write(6,*)'No data are marked as equilibration data.' - endif ! Nequil - endif ! verbose - if(Nlines<2)then - write(6,*)'There are less than two lines of data in '//trim(filename)//'.' - write(6,*)'One cannot obtain error bars with fewer than 2 data points.' - stop - endif - -! Get title. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,'(a)',iostat=ierr)datastring - call check_ierr(ierr) - s=index(datastring,'#') - if(s>0)then - title=datastring(s+1:len_trim(datastring)) - else - write(6,*)'Header line does not have a "#" in front. Stopping.' - stop - endif - title=adjustl(title) - if(verbose)write(6,*)'Title: '//trim(title) - -! Get version number. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,version - call check_ierr(ierr) - call check_hash(temp) - if(verbose)write(6,*)'File version number is '//trim(i2s(version))//'.' - if(version/=1)then - write(6,*)'Version number of '//trim(filename)//' must be 1.' - stop - endif ! version/=1 - -! Get QMC method. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,qmc_method - call check_ierr(ierr) - call check_hash(temp) - qmc_method=adjustl(qmc_method) - if(verbose)write(6,*)'The data were generated using '//trim(qmc_method)//'.' - if(trim(qmc_method)/='VMC'.and.trim(qmc_method)/='DMC')then - write(6,*)'Method in '//trim(filename)//' should be either VMC or DMC.' - stop - endif ! method - if(trim(filename)=='vmc.hist'.and.trim(qmc_method)/='VMC')then - write(6,*)'Warning: you appear to have non-VMC data in a file called & - &vmc.hist.' - write(6,*) - endif - if(trim(filename)=='dmc.hist'.and.trim(qmc_method)/='DMC')then - write(6,*)'Warning: you appear to have non-DMC data in a file called & - &dmc.hist.' - write(6,*) - endif - -! Get interaction-type (interaction). - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,interaction - call check_ierr(ierr) - call check_hash(temp) - coul_ewald=.false. ; coul_mpc=.false. - select case(trim(interaction)) - case('none','coulomb','ewald','mpc','ewald_mpc','mpc_ewald','manual') - continue - case('1') ; interaction='default' - case('2') ; interaction='mpc' - case('3') ; interaction='ewald_mpc' - case('4') ; interaction='mpc_ewald' - case default - write(6,*)'Value of INTERACTION=',trim(interaction),' not recognized. & - &Stopping.' - stop - end select - select case(trim(interaction)) - case('none') ; continue - case('coulomb','ewald','default','manual') ; coul_ewald=.true. - case('mpc') ; coul_mpc=.true. - case('ewald_mpc','mpc_ewald') ; coul_ewald=.true. ; coul_mpc=.true. - end select - if(verbose)write(6,*)'The value of the interaction parameter is ',& - &trim(interaction),'.' - -! Get constant (ion-ion) energy. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,constant_energy - call check_ierr(ierr) - call check_hash(temp) - if(verbose)write(6,*)'Have got constant energy component.' - -! Get total number of electrons. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,netot - call check_ierr(ierr) - call check_hash(temp) - if(verbose)then - if(netot/=1)then - write(6,*)'There are '//trim(i2s(netot))//' particles in the simulation.' - else - write(6,*)'There is 1 particle in the simulation.' - endif - endif ! verbose - if(netot<1)then - write(6,*)'Should be more than one particle!' - stop - endif - -! Get number of atoms per primitive cell. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,nbasis - call check_ierr(ierr) - call check_hash(temp) - if(verbose)then - if(nbasis/=1)then - write(6,*)'The primitive cell contains '//trim(i2s(nbasis))//' atoms.' - else - write(6,*)'The primitive cell contains 1 atom.' - endif - endif ! verbose - if(nbasis<0)then - write(6,*)'There should be at least zero atoms...' - stop - endif - -! Get number of primitive cells. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,npcells - call check_ierr(ierr) - call check_hash(temp) - if(verbose)then - if(npcells/=1)then - write(6,*)'There are '//trim(i2s(npcells))//' primitive cells.' - else - write(6,*)'There is 1 primitive cell.' - endif - endif ! verbose - if(npcells<1)then - write(6,*)'There should be at least one primitive cell.' - stop - endif - -! When forces are present, allocate force array. - if(forces)then - nitot_max_forces=nbasis*npcells - allocate(tag_forces(22,3,nitot_max_forces),stat=ialloc) - if(ialloc/=0)then - write(6,*)'Force array allocation problem.' - stop - endif ! ialloc/=0 - endif ! forces - -! Basis-type keyword. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,atom_basis_type - call check_ierr(ierr) - call check_hash(temp) - select case(trim(atom_basis_type)) - case('0') ; atom_basis_type='none' - case('1') ; atom_basis_type='plane-wave' - case('2') ; atom_basis_type='gaussian' - case('3') ; atom_basis_type='numerical' - case('4') ; atom_basis_type='blip' - case('5') ; atom_basis_type='non_int_he' - case default - continue - end select - if(verbose)write(6,*)'The value of the atom_basis_type parameter is ' & - &//trim(atom_basis_type)//'.' - -! Get periodicity. - read(io,*,iostat=ierr) - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,isper_flag - call check_ierr(ierr) - call check_hash(temp) - if(isper_flag==1)then - isperiodic=.true. ; if(verbose)write(6,*)'The system is periodic.' - select case(trim(interaction)) - case('default','coulomb') ; interaction='ewald' - end select - elseif(isper_flag==0)then - isperiodic=.false. ; if(verbose)write(6,*)'The system is not periodic.' - select case(trim(interaction)) - case('default','ewald') ; interaction='coulomb' - case('mpc','mpc_ewald','ewald_mpc') - write(6,*)'Interaction type should be ''coulomb'' or ''none'' for finite & - &systems. Contradiction in header.' - stop - end select - else - write(6,*)'Periodicity flag must be 0 or 1.' - stop - endif ! periodicity. - -! Get number of data columns. Increase it by 1, since the line-numbers will -! also be read. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - read(io,*,iostat=ierr)temp,no_cols_qmc - call check_ierr(ierr) - call check_hash(temp) - if(verbose)then - if(no_cols_qmc/=1)then - write(6,*)'There are '//trim(i2s(no_cols_qmc))//' columns of data in ' & - &//trim(filename)//'.' - else - write(6,*)'There is 1 column of data in '//trim(filename)//'.' - endif ! Singular/plural - endif ! verbose - if(no_cols_qmc<1)then - write(6,*)'No data to analyse. Stopping.' - stop - endif - no_cols_qmc=no_cols_qmc+1 -! Account for line breaking as the maximum number of items per line is 25 - nbreak=no_cols_qmc/25 - if(modulo(no_cols_qmc,25)>0)nbreak=nbreak+1 - Nlines=Nlines/nbreak - Nequil=Nequil/nbreak - -! Get items in .hist file - tag_step=1 ! Move number - tag_energy=-1 ! Total energy - tag_etotalt=-1 ! Alternative total energy - tag_esqr=-1 ! Square of total energy - tag_popavgsqr=-1 ! Square of population average over total energy - tag_K=-1 ! KEI kinetic-energy estimator - tag_T=-1 ! TI kinetic-energy estimator - tag_fisq=-1 ! FISQ kinetic-energy estimator - tag_Ewald=-1 ! 1/r or Ewald e-e interaction - tag_local=-1 ! Local electron-ion energy - tag_nonlocal=-1 ! Nonlocal electron-ion energy - tag_short=-1 ! Short-range part of MPC - tag_long=-1 ! Long-range part of MPC - tag_cppei=-1 ! Electron-ion CPP term - tag_cppe=-1 ! Electron CPP term - tag_cppee=-1 ! Electron-electron CPP term - tag_masspol=-1 ! Mass-polarization term - tag_future0=-1 ! Future-walking estimator - tag_future1=-1 ! " - tag_future2=-1 ! " - tag_future3=-1 ! " - tag_future4=-1 ! " - tag_future5=-1 ! " - tag_future6=-1 ! " - tag_future7=-1 ! " - tag_future8=-1 ! " - tag_future9=-1 ! " - tag_future10=-1 ! " - tag_massvel=-1 ! Mass-velocity term - tag_darwinen=-1 ! Darwin e-n term - tag_darwinee=-1 ! Darwin e-e term - tag_retard=-1 ! Retardation term. - tag_weight=-1 ! Total weight of configs - tag_nconf=-1 ! Number of configs - tag_eref=-1 ! Reference energy - tag_ebest=-1 ! Best estimate of energy - tag_acc=-1 ! Acceptance ratio - tag_teff=-1 ! Effective time step - tag_dipole1=-1 ! Electric dipole moment - tag_dipole2=-1 ! " " " - tag_dipole3=-1 ! " " " - tag_dipole_sq=-1 ! " " " - tag_contact_den=-1 ! Electron-positron contact density - if(forces)then - tag_forces(1:22,1:3,1:nitot_max_forces)=-1 - nitem_forces=0 ; naxis_forces=0 ; nitot_forces=0 - endif ! forces - - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - do i=2,no_cols_qmc - read(io,*,iostat=ierr)temp,datastring - call check_ierr(ierr) - call check_hash(temp) - datastring=adjustl(datastring) - if(trim(datastring)=='ETOT')then - call check_tag_free(tag_energy) - tag_energy=i - elseif(trim(datastring)=='ETOTALT')then - call check_tag_free(tag_etotalt) - tag_etotalt=i - elseif(trim(datastring)=='ESQR')then - call check_tag_free(tag_esqr) - tag_esqr=i - elseif(trim(datastring)=='POPAVGSQR')then - call check_tag_free(tag_popavgsqr) - tag_popavgsqr=i - elseif(trim(datastring)=='KEI')then - call check_tag_free(tag_K) - tag_K=i - elseif(trim(datastring)=='TI')then - call check_tag_free(tag_T) - tag_T=i - elseif(trim(datastring)=='FISQ')then - call check_tag_free(tag_fisq) - tag_fisq=i - elseif(trim(datastring)=='EWALD')then - call check_tag_free(tag_Ewald) - tag_Ewald=i - elseif(trim(datastring)=='LOCAL')then - call check_tag_free(tag_local) - tag_local=i - elseif(trim(datastring)=='NONLOCAL')then - call check_tag_free(tag_nonlocal) - tag_nonlocal=i - elseif(trim(datastring)=='SHORT')then - call check_tag_free(tag_short) - tag_short=i - elseif(trim(datastring)=='LONG')then - call check_tag_free(tag_long) - tag_long=i - elseif(trim(datastring)=='CPPEI')then - call check_tag_free(tag_cppei) - tag_cppei=i - elseif(trim(datastring)=='CPPE')then - call check_tag_free(tag_cppe) - tag_cppe=i - elseif(trim(datastring)=='CPPEE')then - call check_tag_free(tag_cppee) - tag_cppee=i - elseif(trim(datastring)=='MASSPOL')then - call check_tag_free(tag_masspol) - tag_masspol=i - elseif(trim(datastring(1:3))=='FOR')then - call generate_tag_forces(datastring,i) - elseif(trim(datastring)=='FUTURE0')then - call check_tag_free(tag_future0) - tag_future0=i - elseif(trim(datastring)=='FUTURE1')then - call check_tag_free(tag_future1) - tag_future1=i - elseif(trim(datastring)=='FUTURE2')then - call check_tag_free(tag_future2) - tag_future2=i - elseif(trim(datastring)=='FUTURE3')then - call check_tag_free(tag_future3) - tag_future3=i - elseif(trim(datastring)=='FUTURE4')then - call check_tag_free(tag_future4) - tag_future4=i - elseif(trim(datastring)=='FUTURE5')then - call check_tag_free(tag_future5) - tag_future5=i - elseif(trim(datastring)=='FUTURE6')then - call check_tag_free(tag_future6) - tag_future6=i - elseif(trim(datastring)=='FUTURE7')then - call check_tag_free(tag_future7) - tag_future7=i - elseif(trim(datastring)=='FUTURE8')then - call check_tag_free(tag_future8) - tag_future8=i - elseif(trim(datastring)=='FUTURE9')then - call check_tag_free(tag_future9) - tag_future9=i - elseif(trim(datastring)=='FUTURE10')then - call check_tag_free(tag_future10) - tag_future10=i - elseif(trim(datastring)=='MASSVEL')then - call check_tag_free(tag_massvel) - tag_massvel=i - elseif(trim(datastring)=='DARWINEN')then - call check_tag_free(tag_darwinen) - tag_darwinen=i - elseif(trim(datastring)=='DARWINEE')then - call check_tag_free(tag_darwinee) - tag_darwinee=i - elseif(trim(datastring)=='RETARD')then - call check_tag_free(tag_retard) - tag_retard=i - elseif(trim(datastring)=='WEIGHT')then - call check_tag_free(tag_weight) - tag_weight=i - elseif(trim(datastring)=='NCONF')then - call check_tag_free(tag_nconf) - tag_nconf=i - elseif(trim(datastring)=='EREF')then - call check_tag_free(tag_eref) - tag_eref=i - elseif(trim(datastring)=='EBEST')then - call check_tag_free(tag_ebest) - tag_ebest=i - elseif(trim(datastring)=='ACC')then - call check_tag_free(tag_acc) - tag_acc=i - elseif(trim(datastring)=='TEFF')then - call check_tag_free(tag_teff) - tag_teff=i - elseif(trim(datastring)=='DIPOLE1')then - call check_tag_free(tag_dipole1) - tag_dipole1=i - elseif(trim(datastring)=='DIPOLE2')then - call check_tag_free(tag_dipole2) - tag_dipole2=i - elseif(trim(datastring)=='DIPOLE3')then - call check_tag_free(tag_dipole3) - tag_dipole3=i - elseif(trim(datastring)=='DIPOLESQ')then - call check_tag_free(tag_dipole_sq) - tag_dipole_sq=i - elseif(trim(datastring)=='CONTACT_DEN')then - call check_tag_free(tag_contact_den) - tag_contact_den=i - else - write(6,*)'Column label not recognised.' - write(6,*)'Label is: '//trim(datastring) - stop - endif ! Label - enddo ! i - if(verbose)then - write(6,*)'Have read in column labels.' - write(6,*) - endif ! verbose - -! Warn about missing data, etc. - if(tag_energy<=0)then - write(6,*)'Warning: total energy data are not present!' - write(6,*) - endif - if(tag_K<=0)then - write(6,*)'Warning: kinetic energy (K) data are not present!' - write(6,*) - endif - if((tag_short>0.or.tag_long>0).and..not.coul_mpc)then - write(6,*)'Warning: MPC data are inexplicably present.' - write(6,*) - endif - if(tag_ewald>0.and..not.coul_ewald)then - write(6,*)'Warning: Ewald data are inexplicably present.' - write(6,*) - endif - if(tag_short>0.and.tag_long<=0)then - write(6,*)'Warning: only have short-ranged part of MPC interaction.' - write(6,*) - endif - if(tag_short<=0.and.tag_long>0)then - write(6,*)'Warning: only have long-ranged part of MPC interaction.' - write(6,*) - endif - -! Read final comment line in header. - read(io,*,iostat=ierr)temp - call check_ierr(ierr) - call check_hash(temp) - - END SUBROUTINE read_header - - - SUBROUTINE check_hash(char) -!---------------------------------------------------------------------------! -! This sub is used to check that the 1st char in each header line is a "#". ! -!---------------------------------------------------------------------------! - IMPLICIT NONE - CHARACTER(1),INTENT(in) :: char - if(char/='#')then - write(6,*)'Header line does not have a "#" in front. Stopping.' - stop - endif - END SUBROUTINE check_hash - - - SUBROUTINE check_ierr(ierr,nline) -!------------------------------------------------------! -! Complain if there has been a problem reading a file. ! -!------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: ierr - INTEGER,INTENT(in),OPTIONAL :: nline - if(ierr/=0)then - if(present(nline))then - write(6,*)'Problem reading '//trim(filename)//' at line '//trim(i2s(nline))//'.' - else - write(6,*)'Problem reading '//trim(filename)//'.' - endif - stop - endif - END SUBROUTINE check_ierr - - - SUBROUTINE check_tag_free(tag) -!----------------------------------------------! -! Complain if a tag has already been assigned. ! -!----------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: tag - if(tag/=-1)then - write(6,*)'Tag assigned twice. Two column labels must be the same.' - stop - endif - END SUBROUTINE check_tag_free - - - SUBROUTINE read_data(dmc) -!--------------------------------------------------! -! Read in the raw QMC data from the .hist file. ! -!--------------------------------------------------! - IMPLICIT NONE - LOGICAL,INTENT(in) :: dmc - INTEGER ierr,i,ialloc,nbreak,nleft,in,im - INTEGER,PARAMETER :: io=8 - CHARACTER(640) char640 - -! Open the data file. - open(unit=io,file=trim(filename),status='old',iostat=ierr) - if(ierr/=0)then - write(6,*)'Sorry, cannot open '//trim(filename)//'.' - stop - endif - -! Count the columns and rows of data, establish which data are present, etc. - call read_header(io,dmc) - -! Allocate the data array. - allocate(data_array(Nlines,no_cols_qmc),stat=ialloc) - if(ialloc/=0)then - write(6,*)'Allocation problem (1).' - stop - endif - -! Read in the data. Ignore comments. - i=0 - do - read(io,'(a)',iostat=ierr)char640 - call check_ierr(ierr,i+1) - if(index(char640,'#')==0)then - i=i+1 -! When reading from .hist file, account for maximum number of items -! (=25 per line). - nbreak=no_cols_qmc/25 - nleft=modulo(no_cols_qmc,25) - im=0 - if(nbreak>0)then - do in=1,nbreak - im=in - read(char640,*,iostat=ierr)data_array(i,(im-1)*25+1:im*25) - read(io,'(a)',iostat=ierr)char640 - enddo - endif - if(nleft>0)then - read(char640,*,iostat=ierr)data_array(i,im*25+1:no_cols_qmc) - endif - call check_ierr(ierr,i+1) - if(i>=Nlines)exit - endif ! Line not a comment. - enddo ! i - close(io) - - END SUBROUTINE read_data - - - SUBROUTINE check_data -!------------------------------------------------------------------------! -! This subroutine checks that the raw data in data_array are consistent. ! -! It checks that adding up the energy components gives the total energy, ! -! and that the ion-ion energy in the header is correct. It looks for ! -! Ewald and MPC data and decides which is to be used in the total energy.! -! The number of equilibration steps to be discarded are chosen and the ! -! energy units are selected. ! -!------------------------------------------------------------------------! - IMPLICIT NONE - INTEGER i,ierr,units_choice - DOUBLE PRECISION econst_check,econst_hist,escale,tot_weight - LOGICAL econst_is_const -! Tolerance for checking that total energy is sum of components. - DOUBLE PRECISION,PARAMETER :: tol=1.d-6 - -! Check move numbers - if(tag_step>0)then - do i=1,Nlines - if(nint(data_array(i,tag_step))/=i)then - write(6,*)'WARNING: iteration number behaves oddly at line ' & - &//trim(i2s(i))//' in '//trim(filename)//'.' - write(6,*) - exit - endif ! Problem with move number - enddo ! i - endif ! tag_step>0 - -! Check weights. - if(tag_weight>0)then - tot_weight=0.d0 - do i=1,Nlines - if(data_array(i,tag_weight)<0.d0)then - write(6,*)'Found a negative weight at line '//trim(i2s(i)) & - &//' of '//trim(filename)//'.' - stop - endif ! weight<0 - tot_weight=tot_weight+data_array(i,tag_weight) - enddo ! i - if(tot_weight<=0.d0)then - write(6,*)'Sum of weights is 0. Stopping.' - stop - endif ! total weight=0 - if(.not.use_weights)then - write(6,*)'Weights are present, but will not be used.' - write(6,*) - endif ! weights not to be used. - endif ! weights present. - -! Check that total energy minus KE, e-i pot E, and e-e pot E is a constant: -! the ion-ion energy. - if(tag_energy>0)then - econst_is_const=.true. - do i=1,Nlines - econst_check=data_array(i,tag_energy) - if(tag_K>0)econst_check=econst_check-data_array(i,tag_K) - if(trim(interaction)=='coulomb'.or.trim(interaction)=='ewald'.or.& - &trim(interaction)=='ewald_mpc'.or.trim(interaction)=='manual')then - if(tag_ewald>0)econst_check=econst_check-data_array(i,tag_ewald) - elseif(trim(interaction)=='mpc'.or.trim(interaction)=='mpc_ewald')then - if(tag_short>0)econst_check=econst_check-data_array(i,tag_short) - if(tag_long>0)econst_check=econst_check-data_array(i,tag_long) - endif ! MPC or Ewald present. - if(tag_local>0)econst_check=econst_check-data_array(i,tag_local) - if(tag_nonlocal>0)econst_check=econst_check-data_array(i,tag_nonlocal) - if(tag_cppei>0)econst_check=econst_check-data_array(i,tag_cppei) - if(tag_cppe>0)econst_check=econst_check-data_array(i,tag_cppe) - if(tag_cppee>0)econst_check=econst_check-data_array(i,tag_cppee) - if(i==1)then - econst_hist=econst_check - else - if(abs(econst_hist-econst_check)>tol)then - write(6,*)'Warning: some component of energy is not accounted for!' - write(6,*)'First evaluation of ion-ion energy: ',econst_hist - write(6,*)'Later evaluation of ion-ion energy: ',econst_check - write(6,*) - econst_is_const=.false. - exit - endif ! constanet_energy not constant. - endif ! i=1 - enddo ! i - if(abs(econst_hist)tol)then - write(6,*)'Warning: value of ion-ion energy obtained from raw data & - &differs from the value' - write(6,*)'in the header. Missing constant energy component?' - write(6,*) - endif ! Difference in constant_energy - endif ! Can compare constant_energy values - endif ! Components for check present? - -! Check that FISQ=2*TI-KEI. - if(tag_K>0.and.tag_T>0.and.tag_fisq>0)then - do i=1,Nlines - if(abs(data_array(i,tag_fisq)+data_array(i,tag_K) & - &-2.d0*data_array(i,tag_T))>tol)then - write(6,*)'Warning: problem with kinetic-energy estimators. & - &FISQ /= 2.TI-KEI.' - write(6,*)' KEI at line '//trim(i2s(i))//' : ',data_array(i,tag_K) - write(6,*)' TI at line '//trim(i2s(i))//' : ',data_array(i,tag_T) - write(6,*)'FISQ at line '//trim(i2s(i))//' : ',data_array(i,tag_fisq) - write(6,*) - exit - endif ! Problem with KE estimators. - enddo ! i - endif ! All KE estimators present. - -! Find out how many lines are to be skipped. - if(trim(qmc_method)=='DMC')then - do - write(6,*)'There are '//trim(i2s(Nlines))//' lines of data in total.' - write(6,*)'There are '//trim(i2s(Nequil))//' lines of equilibration data.' - write(6,*)'How many initial lines of data do you wish to discard?' - read(5,*,iostat=ierr)Nskip - if(ierr/=0)Nskip=-1 - if(Nskip<0.or.Nskip>Nlines-2)then - write(6,*)'Number of lines to skip must be between 0 and '// & - &trim(i2s(Nlines-2))//'.' - else - exit - endif ! Problem with Nskip - enddo ! Loop asking for Nskip - if(Nskip=1.and.units_choice<=2)then - exit - else - write(6,*)'Please try again. Choose a number between 1 and 2.' - write(6,*) - endif - enddo ! choice loop - if(units_choice==1)then - escale=1.d0 - e_units='au/particle' - else - escale=htoev - e_units='eV/particle' - endif ! units choice - elseif(isperiodic.and.nbasis>0)then -! Periodic system. - do - write(6,*)'Your data are for a periodic system with atoms.' - write(6,*)'Please select units for your energy data.' - write(6,*)'Choose one of: (1) au per prim cell; (2) eV per & - &prim cell;' - write(6,*)'(3) kcal per prim cell; (4) au per atom; (5) eV per atom; & - &(6) kcal per atom.' - read(5,*,iostat=ierr)units_choice - if(ierr/=0)units_choice=-1 - if(units_choice>=1.and.units_choice<=6)then - exit - else - write(6,*)'Please try again. Choose a number between 1 and 6.' - write(6,*) - endif - enddo ! choice loop - if(units_choice==1)then - escale=1.d0 - e_units='au/prim cell' - elseif(units_choice==2)then - escale=htoev - e_units='eV/prim cell' - elseif(units_choice==3)then - escale=htokcal - e_units='kcal/prim cell' - elseif(units_choice==4)then - escale=1.d0/dble(nbasis) - e_units='au/atom' - elseif(units_choice==5)then - escale=htoev/dble(nbasis) - e_units='eV/atom' - else - escale=htokcal/dble(nbasis) - e_units='kcal/atom' - endif - elseif(.not.isperiodic.and.nbasis>0)then -! Finite system. - do - write(6,*)'Your data are for a finite system with atoms.' - write(6,*)'Please select units for your energy data.' - write(6,*)'Choose one of: (1) au; (2) eV; (3) kcal; (4) au per atom; & - &(5) eV per atom;' - write(6,*)'(6) kcal per atom.' - read(5,*,iostat=ierr)units_choice - if(ierr/=0)units_choice=-1 - if(units_choice>=1.and.units_choice<=6)then - exit - else - write(6,*)'Please try again. Choose a number between 1 and 6.' - write(6,*) - endif - enddo ! choice loop - if(units_choice==1)then - escale=1.d0 - e_units='au' - elseif(units_choice==2)then - escale=htoev - e_units='eV' - elseif(units_choice==3)then - escale=htokcal - e_units='kcal' - elseif(units_choice==4)then - escale=1.d0/dble(nbasis) - e_units='au/atom' - elseif(units_choice==5)then - escale=htoev/dble(nbasis) - e_units='eV/atom' - else - escale=htokcal/dble(nbasis) - e_units='kcal/atom' - endif - else -! Default. - do - write(6,*)'Please select units for your energy data.' - write(6,*)'Choose one of: (1) au; (2) eV; (3) kcal.' - read(5,*,iostat=ierr)units_choice - if(ierr/=0)units_choice=-1 - if(units_choice>=1.and.units_choice<=3)then - exit - else - write(6,*)'Please try again. Choose a number between 1 and 3.' - write(6,*) - endif - enddo ! choice loop - if(units_choice==1)then - escale=1.d0 - e_units='au' - elseif(units_choice==2)then - escale=htoev - e_units='eV' - else - escale=htokcal - e_units='kcal' - endif - write(6,*)'For finite systems, energies are quoted for the whole system.' - write(6,*)'For real crystals, energies are quoted per primitive cell.' - write(6,*)'For electron(-hole) systems, energies are quoted per particle.' - endif - write(6,*) - -! Rescale energy data. (Easier just to rescale the raw data than to rescale -! each result quoted.) - if(escale/=1.d0)then - if(tag_energy>0)data_array(:,tag_energy)=data_array(:,tag_energy)*escale - if(tag_etotalt>0)data_array(:,tag_etotalt)=data_array(:,tag_etotalt)*escale - if(tag_esqr>0)data_array(:,tag_esqr)=data_array(:,tag_esqr)*escale*escale - if(tag_popavgsqr>0)data_array(:,tag_popavgsqr)=data_array(:,tag_popavgsqr) & - &*escale*escale - if(tag_K>0)data_array(:,tag_K)=data_array(:,tag_K)*escale - if(tag_T>0)data_array(:,tag_T)=data_array(:,tag_T)*escale - if(tag_fisq>0)data_array(:,tag_fisq)=data_array(:,tag_fisq)*escale - if(tag_Ewald>0)data_array(:,tag_Ewald)=data_array(:,tag_Ewald)*escale - if(tag_local>0)data_array(:,tag_local)=data_array(:,tag_local)*escale - if(tag_nonlocal>0)data_array(:,tag_nonlocal)=data_array(:,tag_nonlocal) & - &*escale - if(tag_short>0)data_array(:,tag_short)=data_array(:,tag_short)*escale - if(tag_long>0)data_array(:,tag_long)=data_array(:,tag_long)*escale - if(tag_cppei>0)data_array(:,tag_cppei)=data_array(:,tag_cppei)*escale - if(tag_cppe>0)data_array(:,tag_cppe)=data_array(:,tag_cppe)*escale - if(tag_cppee>0)data_array(:,tag_cppee)=data_array(:,tag_cppee)*escale - if(tag_masspol>0)data_array(:,tag_masspol)=data_array(:,tag_masspol)*escale - if(forces)then - do iion=1,nitot_forces - do iaxis=1,3 - do item=1,22 - if(tag_forces(item,iaxis,iion)/=-1)then - data_array(:,tag_forces(item,iaxis,iion))=& - &data_array(:,tag_forces(item,iaxis,iion))*escale - endif - enddo - enddo - enddo - endif ! if forces - if(tag_future0>0)data_array(:,tag_future0)=data_array(:,tag_future0)*escale - if(tag_future1>0)data_array(:,tag_future1)=data_array(:,tag_future1)*escale - if(tag_future2>0)data_array(:,tag_future2)=data_array(:,tag_future2)*escale - if(tag_future3>0)data_array(:,tag_future3)=data_array(:,tag_future3)*escale - if(tag_future4>0)data_array(:,tag_future4)=data_array(:,tag_future4)*escale - if(tag_future5>0)data_array(:,tag_future5)=data_array(:,tag_future5)*escale - if(tag_future6>0)data_array(:,tag_future6)=data_array(:,tag_future6)*escale - if(tag_future7>0)data_array(:,tag_future7)=data_array(:,tag_future7)*escale - if(tag_future8>0)data_array(:,tag_future8)=data_array(:,tag_future8)*escale - if(tag_future9>0)data_array(:,tag_future9)=data_array(:,tag_future9)*escale - if(tag_future10>0)data_array(:,tag_future10)=data_array(:,tag_future10)*& - &escale - if(tag_massvel>0)data_array(:,tag_massvel)=data_array(:,tag_massvel)*escale - if(tag_darwinen>0)data_array(:,tag_darwinen)=data_array(:,tag_darwinen) & - &*escale - if(tag_darwinee>0)data_array(:,tag_darwinee)=data_array(:,tag_darwinee) & - &*escale - if(tag_retard>0)data_array(:,tag_retard)=data_array(:,tag_retard)*escale - if(tag_eref>0)data_array(:,tag_eref)=data_array(:,tag_eref)*escale - if(tag_ebest>0)data_array(:,tag_ebest)=data_array(:,tag_ebest)*escale - constant_energy=constant_energy*escale - endif ! Data needs rescaling. - - END SUBROUTINE check_data - - - SUBROUTINE compute_stats -!--------------------------------------------------------------------------! -! In this subroutine, the various columns of data are subjected to various ! -! statistical analyses. ! -!--------------------------------------------------------------------------! - IMPLICIT NONE - INTEGER ierr,block_length,Nstudy,startline,nthird,nthirdstart, & - &ialloc,nthirdstop,i - DOUBLE PRECISION av,av_energy,std_err,std_err_energy,delta_std_err,var, & - &max_val,min_val,skew,kurt,corr_tau,corr_tau_err,sqrt_tau,err_sqrt_tau, & - &raw_var,raw_var_err,pop_var,pop_var_err - DOUBLE PRECISION,ALLOCATABLE :: temp_data(:) - - Nstudy=Nlines-Nskip - startline=Nskip+1 - -! Write out some information about the more important DMC simulation params. -! Do this first, so that important data appears at end of output. - if(tag_nconf>0)then - - write(6,*)'ANALYSIS OF CONFIGURATION POPULATION' - write(6,*)'====================================' - call compute_stats_unweighted(.false.,Nstudy,data_array(startline:Nlines, & - &tag_nconf),av,var,skew,kurt,max_val,min_val) - write(6,*)'Minimum population : ',min_val - write(6,*)' Mean population : ',av - write(6,*)'Maximum population : ',max_val - write(6,*)' Std error : ',sqrt(var/dble(Nstudy)) - if(av-min_val>0.25d0*av.or.max_val-av>0.25d0*av)write(6,*) & - &'Warning: Population fluctuated by more than 25% of mean.' - write(6,*) - - endif ! Config population data present. - - if(tag_acc>0)then - - write(6,*)'ANALYSIS OF ACCEPTANCE RATIO' - write(6,*)'============================' - call compute_stats_unweighted(.false.,Nstudy,data_array(startline:Nlines, & - &tag_acc),av,var,skew,kurt,max_val,min_val) - write(6,*)'Minimum acceptance ratio : ',min_val - write(6,*)' Mean acceptance ratio : ',av - write(6,*)'Maximum acceptance ratio : ',max_val - write(6,*)' Std error : ',sqrt(var/dble(Nstudy)) - write(6,*) - - endif ! Acceptance-ratio data present. - - if(tag_teff>0)then - - write(6,*)'ANALYSIS OF EFFECTIVE TIME STEP' - write(6,*)'===============================' - call compute_stats_unweighted(.false.,Nstudy,data_array(startline:Nlines, & - &tag_teff),av,var,skew,kurt,max_val,min_val) - write(6,*)'Minimum time step (au) : ',min_val - write(6,*)' Mean time step (au) : ',av - write(6,*)'Maximum time step (au) : ',max_val - write(6,*)' Std error (au) : ',sqrt(var/dble(Nstudy)) - write(6,*) - - endif ! Effective time step data present. - - if(tag_energy>0)then - -! Compute lots of information about the total energy data. - write(6,*)'ANALYSIS OF TOTAL-ENERGY DATA' - write(6,*)'=============================' - call compute_stats_unweighted(.true.,Nstudy,data_array(startline:Nlines, & - &tag_energy),av,var,skew,kurt,max_val,min_val) - write(6,*)'Minimum energy (',trim(e_units),') : ',min_val - write(6,*)' Mean energy (',trim(e_units),') : ',av - write(6,*)'Maximum energy (',trim(e_units),') : ',max_val - write(6,*)' Variance (',trim(e_units),') : ',var - write(6,*)' Std error (',trim(e_units),') : ',sqrt(var/dble(Nstudy)) - write(6,*)repeat(' ',len_trim(e_units))//' Skewness : ',skew - write(6,*)repeat(' ',len_trim(e_units))//'Normal sk. fluct. : ',& - &sqrt(6.d0/dble(Nstudy)) - write(6,*)repeat(' ',len_trim(e_units))//' Kurtosis : ',kurt - write(6,*)repeat(' ',len_trim(e_units))//'Normal ku. fluct. : ',& - &sqrt(24.d0/dble(Nstudy)) - write(6,*)'(NB, the var of the energy data is not an estimate & - &of the actual var.)' - write(6,*) -! Analyse total energy by thirds if there is enough data. - if(Nstudy>=6)then - write(6,*)'ANALYSIS OF TOTAL-ENERGY DATA BY THIRDS' - write(6,*)'=======================================' - write(6,*)'(Energy data in units of '//trim(e_units)//'.)' - nthird=Nstudy/3 - write(6,*)' Data range Av energy Std error Maximum & - &Minimum' - do i=1,3 - nthirdstart=startline+(i-1)*nthird - nthirdstop=nthirdstart+nthird-1 - call compute_stats_unweighted(.false.,nthird, & - &data_array(nthirdstart:nthirdstop,tag_energy),av,var,skew,kurt, & - &max_val,min_val) - write(6,'(" ",a16,4(" ",es12.5))')trim(i2s(nthirdstart)) & - &//'->'//trim(i2s(nthirdstop)),av,sqrt(var/dble(nthird)),max_val,min_val - enddo ! i - else - write(6,*)'Not enough data to analyse by thirds: need at least 6 points.' - endif ! Enough data? - write(6,*) - write(6,*)'CORRELATION-TIME ANALYSIS OF TOTAL-ENERGY DATA' - write(6,*)'==============================================' - call correlation_time(Nstudy,data_array(startline:Nlines,tag_energy), & - &corr_tau,corr_tau_err,av,var) - if(corr_tau/=-1.d0)then - write(6,*)' Correlation time (steps) : ',corr_tau - write(6,*)' Error in correlation time (steps) : ',corr_tau_err - write(6,*) - if(corr_tau>0.d0)then - sqrt_tau=sqrt(corr_tau) - err_sqrt_tau=corr_tau_err/(2*sqrt_tau) - write(6,*)' Error-bar factor : ',sqrt_tau - write(6,*)' Error in error-bar factor : ',err_sqrt_tau - write(6,*) - if(tag_weight>0.and.use_weights)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_energy), & - &data_array(startline:Nlines,tag_weight),1,av,std_err,delta_std_err) - else - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_energy), & - &1,av,std_err,delta_std_err) - endif - write(6,*)' Mean energy (',trim(e_units),') : ',av - write(6,*)' Correlation-corrected error (',trim(e_units),') : ',& - &std_err*sqrt(corr_tau) - write(6,*)' Error in error (',trim(e_units),') : ',& - &sqrt((std_err*err_sqrt_tau)**2+(sqrt_tau*delta_std_err)**2) - else - write(6,*)'The correlation time appears to be negative.' - endif - else - write(6,*)'The correlation time could not be computed.' - endif ! corr_tau calculated. - write(6,*) - - write(6,*)'REBLOCKING ANALYSIS OF TOTAL-ENERGY DATA' - write(6,*)'========================================' -! Print out reblocking analysis of energy - write(6,*)'(Energy data in units of '//trim(e_units)//'.)' - if(tag_weight>0.and.use_weights)then - call reblock_analysis(Nstudy,data_array(startline:Nlines,tag_energy), & - &data_array(startline:Nlines,tag_weight)) - else - call reblock_analysis(Nstudy,data_array(startline:Nlines,tag_energy)) - endif ! weights - - endif ! energy data available - - do - write(6,*)'Please choose a block length for reblocking all energy & - &components.' - read(5,*,iostat=ierr)block_length - if(ierr/=0)block_length=-1 - if(block_length>=1.and.block_length<=Nstudy/2)then - exit - else - write(6,*)'Please try again. Block length should be between 1 and ' & - &//trim(i2s(Nstudy/2))//'.' - endif - enddo ! get block length - write(6,*)'Chosen block length: '//trim(i2s(block_length))//'.' - write(6,*) - -! Write out the energy components with reblocked error bars. - - write(6,*)'ENERGY COMPONENTS WITH REBLOCKED ERROR BARS' - write(6,*)'===========================================' - -5 format(" ",a30,2(" ",es22.14)) -10 format(" ",a30," ",es22.14) -15 format(32x,a23,a) - write(6,15)' Mean ('//trim(e_units)//') ',' Err (' & - &//trim(e_units)//')' - - allocate(temp_data(startline:Nlines),stat=ialloc) - if(ialloc/=0)then - write(6,*)'Allocation problem.' - stop - endif - - if(tag_weight>0.and.use_weights)then - - if(tag_energy>0)then - - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_energy), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - av_energy=av ; std_err_energy=std_err - if(trim(interaction)=='mpc'.or.trim(interaction)=='mpc_ewald')then - write(6,5)'Total energy (using MPC) :',av,std_err - if(trim(interaction)=='mpc_ewald'.and.tag_ewald>0.and.tag_short>0.and.& - &tag_long>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_energy) & - &-data_array(startline:Nlines,tag_short)-data_array(startline:Nlines, & - &tag_long)+data_array(startline:Nlines,tag_ewald) - call reblock_weighted(Nstudy,temp_data(startline:Nlines), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Total energy (using Ewald) :',av,std_err - endif ! ewald present - else - if(isperiodic)then - write(6,5)'Total energy (using Ewald) :',av,std_err - else - write(6,5)'Total energy :',av,std_err - endif ! periodic - if(trim(interaction)=='ewald_mpc'.and.tag_ewald>0.and.tag_short>0.and.& - &tag_long>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_energy) & - &+data_array(startline:Nlines,tag_short)+data_array(startline:Nlines, & - &tag_long)-data_array(startline:Nlines,tag_ewald) - call reblock_weighted(Nstudy,temp_data(startline:Nlines), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Total energy (using MPC) :',av,std_err - endif ! MPC present - endif ! use_mpc_energy - - if(tag_esqr>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_esqr), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Total energy squared:',av,std_err - raw_var=av-av_energy*av_energy - raw_var_err=sqrt(std_err**2+(av_energy*std_err_energy)**2) - write(6,5)'Variance of total energy:',raw_var,raw_var_err - endif - - if(tag_popavgsqr>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_popavgsqr), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Population avg of total energy squared:',av,std_err - pop_var=av-av_energy*av_energy - pop_var_err=sqrt(std_err**2+(av_energy*std_err_energy)**2) - write(6,5)'Variance of population avg of energy:',pop_var,pop_var_err - endif - - if(tag_esqr>0.and.tag_popavgsqr>0)then - write(6,5)'Effective population size:',raw_var/pop_var,& - &sqrt((raw_var_err/pop_var)**2+(pop_var_err*raw_var/pop_var**2)**2) - endif - - if(tag_masspol>0.and.tag_massvel>0.and.tag_darwinen>0.and.tag_darwinee>0 & - &.and.tag_retard>0)then -! At present, only have relativistic data for atoms. - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_energy) & - &+data_array(startline:Nlines,tag_masspol) & - &+data_array(startline:Nlines,tag_massvel) & - &+data_array(startline:Nlines,tag_darwinen) & - &+data_array(startline:Nlines,tag_darwinee) & - &+data_array(startline:Nlines,tag_retard) - call reblock_weighted(Nstudy,temp_data(startline:Nlines), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Total energy (inc rel) :',av,std_err - endif ! rel_present - - if(tag_K>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_energy) & - &-data_array(startline:Nlines,tag_K) - call reblock_weighted(Nstudy,temp_data(startline:Nlines), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - if(trim(interaction)=='mpc'.or.trim(interaction)=='mpc_ewald')then - write(6,5)'Total pot energy (using MPC) :',av,std_err - else - if(isperiodic)then - write(6,5)'Tot pot energy (using Ewald) :',av,std_err - else - write(6,5)'Total potential energy :',av,std_err - endif ! periodic - endif ! use_mpc_energy - endif ! K present - - endif ! energy present. - - if(tag_K>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_K), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Kinetic energy (K) :',av,std_err - endif ! K present - if(tag_T>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_T), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Kinetic energy (T) :',av,std_err - endif ! T present - if(tag_fisq>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_fisq), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Kinetic energy (FISQ) :',av,std_err - elseif(tag_K>0.and.tag_T>0)then - temp_data(startline:Nlines)=2.d0*data_array(startline:Nlines,tag_T) & - &-data_array(startline:Nlines,tag_K) - call reblock_weighted(Nstudy,temp_data(startline:Nlines), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Kinetic energy (FISQ) :',av,std_err - endif ! FISQ present. - if(tag_ewald>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_ewald), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - if(isperiodic)then - write(6,5)'Ewald interaction :',av,std_err - else - write(6,5)'Coulomb interaction :',av,std_err - endif ! periodic - endif ! Ewald present. - if(tag_local>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_local), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Local e-i energy :',av,std_err - endif ! local present - if(tag_nonlocal>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines, & - &tag_nonlocal),data_array(startline:Nlines,tag_weight),block_length, & - &av,std_err,delta_std_err) - write(6,5)'Nonlocal e-i energy :',av,std_err - endif ! nonlocal e-i pot present - - else - - if(tag_energy>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_energy), & - &block_length,av,std_err,delta_std_err) - av_energy=av ; std_err_energy=std_err - if(trim(interaction)=='mpc'.or.trim(interaction)=='mpc_ewald')then - write(6,5)'Total energy (using MPC) :',av,std_err - if(trim(interaction)=='mpc_ewald'.and.tag_ewald>0.and.tag_long>0.and.& - &tag_short>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_energy) & - &-data_array(startline:Nlines,tag_short)-data_array(startline:Nlines, & - &tag_long)+data_array(startline:Nlines,tag_ewald) - call reblock_unweighted(Nstudy,temp_data(startline:Nlines), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Total energy (using Ewald) :',av,std_err - endif ! Ewald present. - else - if(isperiodic)then - write(6,5)'Total energy (using Ewald) :',av,std_err - else - write(6,5)'Total energy :',av,std_err - endif ! periodic - if(trim(interaction)=='ewald_mpc'.and.tag_ewald>0.and.tag_long>0.and.& - &tag_short>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_energy) & - &+data_array(startline:Nlines,tag_short)+data_array(startline:Nlines, & - &tag_long)-data_array(startline:Nlines,tag_ewald) - call reblock_unweighted(Nstudy,temp_data(startline:Nlines), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Total energy (using MPC) :',av,std_err - endif ! MPC present. - endif ! use_mpc_energy - if(tag_esqr>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_esqr), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Total energy squared:',av,std_err - write(6,10)'Raw variance of total energy:',& - &(av-av_energy*av_energy)*dble(nparts_per_simcell) - endif - if(tag_masspol>0.and.tag_massvel>0.and.tag_darwinen>0.and.tag_darwinee>0 & - &.and.tag_retard>0)then -! At present, only have relativistic data for atoms. - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_energy) & - &+data_array(startline:Nlines,tag_masspol) & - &+data_array(startline:Nlines,tag_massvel) & - &+data_array(startline:Nlines,tag_darwinen) & - &+data_array(startline:Nlines,tag_darwinee) & - &+data_array(startline:Nlines,tag_retard) - call reblock_unweighted(Nstudy,temp_data(startline:Nlines),block_length, & - &av,std_err,delta_std_err) - write(6,5)'Total energy (inc rel) :',av,std_err - endif ! rel_present - if(tag_K>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_energy) & - &-data_array(startline:Nlines,tag_K) - call reblock_unweighted(Nstudy,temp_data(startline:Nlines),block_length, & - &av,std_err,delta_std_err) - if(trim(interaction)=='mpc'.or.trim(interaction)=='mpc_ewald')then - write(6,5)'Total pot energy (using MPC) :',av,std_err - else - if(isperiodic)then - write(6,5)'Tot pot energy (using Ewald) :',av,std_err - else - write(6,5)'Total potential energy :',av,std_err - endif ! periodic - endif ! use_mpc_energy - endif ! K present - endif ! Energy data present. - if(tag_K>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_K), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Kinetic energy (K) :',av,std_err - endif ! K present. - if(tag_T>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_T), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Kinetic energy (T) :',av,std_err - endif ! T present - if(tag_fisq>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_fisq), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Kinetic energy (FISQ) :',av,std_err - elseif(tag_K>0.and.tag_T>0)then - temp_data(startline:Nlines)=2.d0*data_array(startline:Nlines,tag_T) & - &-data_array(startline:Nlines,tag_K) - call reblock_unweighted(Nstudy,temp_data(startline:Nlines),block_length, & - &av,std_err,delta_std_err) - write(6,5)'Kinetic energy (FISQ) :',av,std_err - endif ! K & T present. - if(tag_ewald>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_ewald), & - &block_length,av,std_err,delta_std_err) - if(isperiodic)then - write(6,5)'Ewald interaction :',av,std_err - else - write(6,5)'Coulomb interaction :',av,std_err - endif ! periodic - endif ! Ewald present. - if(tag_local>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_local), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Local e-i energy :',av,std_err - endif ! local present - if(tag_nonlocal>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines, & - &tag_nonlocal),block_length,av,std_err,delta_std_err) - write(6,5)'Nonlocal e-i energy :',av,std_err - endif ! nonlocal e-i pot present - - endif ! weighted - - if(tag_weight>0.and.use_weights)then - if(tag_short>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_short), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Short-range MPC energy :',av,std_err - endif ! short present - if(tag_long>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_long), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Long-range MPC energy :',av,std_err - endif ! long present. - if(tag_short>0.and.tag_long>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_short)+ & - &data_array(startline:Nlines,tag_long) - call reblock_weighted(Nstudy,temp_data(startline:Nlines), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Total MPC energy :',av,std_err - endif ! MPC data present - if(tag_cppei>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_cppei), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'CPP energy (e-i) :',av,std_err - endif ! CPPEI present - if(tag_cppe>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_cppe), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'CPP energy (e) :',av,std_err - endif ! CPPE present - if(tag_cppee>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_cppee), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'CPP energy (e-e) :',av,std_err - endif ! CPP data present - if(tag_masspol>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_masspol), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Mass-polarization energy :',av,std_err - endif ! masspol present - if(tag_massvel>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_massvel), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Mass-velocity energy :',av,std_err - endif ! massvel present. - if(tag_darwinen>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_darwinen), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Electron-nucleus Darwin :',av,std_err - endif ! darawin e-n present. - if(tag_darwinee>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_darwinee), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Electron-electron Darwin :',av,std_err - endif ! darwin e-e present. - if(tag_retard>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_retard), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Retardation term :',av,std_err - endif ! Duh. - if(tag_masspol>0.and.tag_massvel>0.and.tag_darwinen>0.and.tag_darwinee>0 & - &.and.tag_retard>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_masspol) & - &+data_array(startline:Nlines,tag_massvel) & - &+data_array(startline:Nlines,tag_darwinen) & - &+data_array(startline:Nlines,tag_darwinee) & - &+data_array(startline:Nlines,tag_retard) - call reblock_weighted(Nstudy,temp_data(startline:Nlines), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Total rel correction :',av,std_err - endif ! Relativistic data present - - else - - if(tag_short>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_short), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Short-range MPC energy :',av,std_err - endif ! short present. - if(tag_long>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_long), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Long-range MPC energy :',av,std_err - endif ! long present. - if(tag_short>0.and.tag_long>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_short) & - &+data_array(startline:Nlines,tag_long) - call reblock_unweighted(Nstudy,temp_data(startline:Nlines),block_length, & - &av,std_err,delta_std_err) - write(6,5)'Total MPC energy :',av,std_err - endif ! MPC data present - if(tag_cppei>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_cppei), & - &block_length,av,std_err,delta_std_err) - write(6,5)'CPP energy (e-i) :',av,std_err - endif ! CPPEI present - if(tag_cppe>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_cppe), & - &block_length,av,std_err,delta_std_err) - write(6,5)'CPP energy (e) :',av,std_err - endif ! CPPE present. - if(tag_cppee>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_cppee), & - &block_length,av,std_err,delta_std_err) - write(6,5)'CPP energy (e-e) :',av,std_err - endif ! CPPEE data present - if(tag_masspol>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_masspol), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Mass-polarization energy :',av,std_err - endif ! masspol present. - if(tag_massvel>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_massvel), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Mass-velocity energy :',av,std_err - endif ! massvel present. - if(tag_darwinen>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_darwinen), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Electron-nucleus Darwin :',av,std_err - endif ! darwinen present. - if(tag_darwinee>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_darwinee), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Electron-electron Darwin :',av,std_err - endif ! darwinee present. - if(tag_retard>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_retard), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Retardation term :',av,std_err - endif ! Duh. - if(tag_masspol>0.and.tag_massvel>0.and.tag_darwinen>0.and.tag_darwinee>0 & - &.and.tag_retard>0)then - temp_data(startline:Nlines)=data_array(startline:Nlines,tag_masspol) & - &+data_array(startline:Nlines,tag_massvel) & - &+data_array(startline:Nlines,tag_darwinen) & - &+data_array(startline:Nlines,tag_darwinee) & - &+data_array(startline:Nlines,tag_retard) - call reblock_unweighted(Nstudy,temp_data(startline:Nlines), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Total rel correction :',av,std_err - endif ! Relativistic data present - - endif ! weighted. - - deallocate(temp_data) - - if(constant_energy/=0.d0)write(6,10)'Constant energy :',constant_energy - write(6,*) - -! Calculate and write out forces - if(forces)call construct_write_forces(startline,Nstudy,block_length) - -! Write out future-walking estimates - if(tag_future1>0)then - write(6,*)'FUTURE-WALKING ESTIMATES WITH REBLOCKED ERROR BAR' - write(6,*)'=================================================' - write(6,*)'Future-walking estimates of the observable pureitems(1) & - &from the dmc_main' - write(6,*)'routine. Temporarily, pureitems(1) is the Hellmann-Feynman & - &forces in the' - write(6,*)'x-direction of the 1st atom as ordered in the gwfn.data & - &file. To estimate' - write(6,*)'a different observable, alter the assignment after line:' - write(6,*)"'Change next line when future-walking estimates are required'" - write(6,*)' Mean (au) Err (au)' - - if(tag_weight>0.and.use_weights)then - if(tag_future0>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future0), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (1st) :',av,std_err - endif - if(tag_future1>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future1), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (2nd) :',av,std_err - endif - if(tag_future2>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future2), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (3rd) :',av,std_err - endif - if(tag_future3>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future3), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (4th) :',av,std_err - endif - if(tag_future4>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future4), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (5th) :',av,std_err - endif - if(tag_future5>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future5), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (6th) :',av,std_err - endif - if(tag_future6>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future6), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (7th) :',av,std_err - endif - if(tag_future7>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future7), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (8th) :',av,std_err - endif - if(tag_future8>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future8), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (9th) :',av,std_err - endif - if(tag_future9>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future9), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (10th) :',av,std_err - endif - if(tag_future10>0) then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_future10), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'FW Estimator (11th) :',av,std_err - endif - write(6,*) - write(6,*)'The 1st estimator corresponds to future-walking projection & - &time T=0 1/Ha,' - write(6,*)'the 2nd to T=0.5 1/Ha, the 3rd to T=1 1/Ha..., and the & - &11th to T=10 1/Ha.' - write(6,*) - - endif ! tag_weight - - endif ! future - - if(tag_dipole1>0.or.tag_dipole2>0.or.tag_dipole3>0.or.tag_dipole_sq>0)then - write(6,*)'ELECTRIC DIPOLE MOMENT WITH REBLOCKED ERROR BARS' - write(6,*)'================================================' - write(6,15)' Mean (au) ',' Err (au)' - - if(tag_weight>0.and.use_weights)then - if(tag_dipole1>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_dipole1), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Dipole moment (x cpt) :',av,std_err - endif ! tag_dipole1 - if(tag_dipole2>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_dipole2), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Dipole moment (y cpt) :',av,std_err - endif ! tag_dipole2 - if(tag_dipole3>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_dipole3), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Dipole moment (z cpt) :',av,std_err - endif ! tag_dipole3 - if(tag_dipole_sq>0)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_dipole_sq), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Dipole moment squared :',av,std_err - endif ! tag_dipole_sq - else - if(tag_dipole1>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_dipole1), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Dipole moment (x cpt) :',av,std_err - endif ! tag_dipole1 - if(tag_dipole2>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_dipole2), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Dipole moment (y cpt) :',av,std_err - endif ! tag_dipole2 - if(tag_dipole3>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines,tag_dipole3), & - &block_length,av,std_err,delta_std_err) - write(6,5)'Dipole moment (z cpt) :',av,std_err - endif ! tag_dipole3 - if(tag_dipole_sq>0)then - call reblock_unweighted(Nstudy,data_array(startline:Nlines, & - &tag_dipole_sq),block_length,av,std_err,delta_std_err) - write(6,5)'Dipole moment squared :',av,std_err - endif ! tag_dipole_sq - endif ! weighted. - write(6,*) - endif ! Dipole moment - - if(tag_contact_den>0)then - write(6,*)'CONTACT DENSITY' - write(6,*)'===============' - write(6,15)' Mean (au) ',' Err (au)' - if(tag_weight>0.and.use_weights)then - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_contact_den), & - &data_array(startline:Nlines,tag_weight),block_length,av,std_err, & - &delta_std_err) - write(6,5)'Elec-pos contact density :',av,std_err - else - call reblock_unweighted(Nstudy,data_array(startline:Nlines, & - &tag_contact_den),block_length,av,std_err,delta_std_err) - write(6,5)'Elec-pos contact density :',av,std_err - endif ! weighted. - write(6,*) - endif ! Contact density - - END SUBROUTINE compute_stats - - - SUBROUTINE reblock_analysis(no_pts,data_array,weight_array) -!--------------------------------------------------------------! -! Compute the weighted average of the data, and calculate the ! -! error bar as a function of reblocking transformation number. ! -!--------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: no_pts - DOUBLE PRECISION,INTENT(in) :: data_array(no_pts) - DOUBLE PRECISION,INTENT(in),OPTIONAL :: weight_array(no_pts) - INTEGER no_rtns,rtn,block_length,ierr - DOUBLE PRECISION av,std_err,delta_std_err - - open(unit=10,file='reblock.plot',status='replace',iostat=ierr) - if(ierr/=0)then - write(6,*)'Problem opening reblock.plot.' - stop - endif - -! Number of reblocking transformations - no_rtns=floor(log(dble(no_pts))/log(2.d0)) - -! Initial block length - block_length=1 - -! Write out results of reblocking analysis - write(6,*)' RTN Blk leng Std error in mean Error in std error' - do rtn=0,no_rtns-1 - if(present(weight_array))then - call reblock_weighted(no_pts,data_array,weight_array,block_length, & - &av,std_err,delta_std_err) - else - call reblock_unweighted(no_pts,data_array,block_length,av,std_err, & - &delta_std_err) - endif ! weights present - write(6,'(" ",i4," ",i10," ",es23.15," ",es23.15)')rtn,block_length, & - &std_err,delta_std_err - write(10,*)rtn,std_err,delta_std_err - block_length=2*block_length - enddo ! rtn - write(6,*) - - write(6,*)'Reblocked error bar against reblocking transformation number & - &(RTN) has been' - write(6,*)'written to reblock.plot. Please use "plot_reblock" to view & - &these data.' - write(6,*) - close(10) - - END SUBROUTINE reblock_analysis - - - SUBROUTINE reblock_forces_analysis(no_pts,data_array,plotname,weight_array) -!--------------------------------------------------------------! -! Compute the weighted average of the data, and calculate the ! -! error bar as a function of reblocking transformation number. ! -! This routine is an extension of routine reblock_analysis ! -! to allow specifying the name of the file to be written out. ! -! ! -! AB 11.2007 ! -!--------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: no_pts - DOUBLE PRECISION,INTENT(in) :: data_array(no_pts) - DOUBLE PRECISION,INTENT(in),OPTIONAL :: weight_array(no_pts) - CHARACTER(20),INTENT(in) :: plotname - INTEGER no_rtns,rtn,block_length,ierr - DOUBLE PRECISION av,std_err,delta_std_err - - open(unit=10,file=plotname ,status='replace',iostat=ierr) - if(ierr/=0)then - write(6,*)'Problem opening plotname.' - stop - endif - -! Number of reblocking transformations - no_rtns=floor(log(dble(no_pts))/log(2.d0)) - -! Initial block length - block_length=1 - -! Write out results of reblocking analysis - write(6,*)' RTN Blk leng Std error in mean Error in std error' - do rtn=0,no_rtns-1 - if(present(weight_array))then - call reblock_weighted(no_pts,data_array,weight_array,block_length, & - &av,std_err,delta_std_err) - else - call reblock_unweighted(no_pts,data_array,block_length,av,std_err, & - &delta_std_err) - endif ! weights present - write(6,'(" ",i4," ",i10," ",es23.15," ",es23.15)')rtn,block_length, & - &std_err,delta_std_err - write(10,*)rtn,std_err,delta_std_err - block_length=2*block_length - enddo ! rtn - write(6,*) - - close(10) - - END SUBROUTINE reblock_forces_analysis - - - CHARACTER(12) FUNCTION i2s(n) -!-----------------------------------------------------------------------! -! I2S ! -! === ! -! Convert integers to left justified strings that can be printed in the ! -! middle of a sentence without introducing large amounts of white space.! -! ! -! Calling routine is intended to include something like: ! -! USE utilities ! -! INTEGER i ! -! i=12 ! -! write(6,*)'Integer number ',trim(i2s(i)),' with words at the end.' ! -!-----------------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: n - INTEGER i,j - INTEGER,PARAMETER :: ichar0=ichar('0') - i2s='' - i=abs(n) - do j=len(i2s),1,-1 - i2s(j:j)=achar(ichar0+mod(i,10)) - i=i/10 ; if(i==0)exit - enddo ! j - if(n<0)then - i2s='-'//adjustl(i2s) - else - i2s=adjustl(i2s) - endif ! n<0 - END FUNCTION i2s - - - SUBROUTINE generate_tag_forces(datastring,i) -!------------------------------------------------------------------! -! This subroutine generates tags for forces from tags read in from ! -! the .hist file. It also determines the number of atoms and axes ! -! for which forces data are available. ! -! ! -! AB 11.2007 ! -!------------------------------------------------------------------! - IMPLICIT NONE - CHARACTER(72),INTENT(in) :: datastring - INTEGER,INTENT(in) :: i - CHARACTER(1) :: axis(3)=(/'X','Y','Z'/) - CHARACTER(1) :: atem(22)=(/'A','B','C','D','E','F','G','H','I','J','K',& - &'L','M','N','O','P','Q','R','S','T','U','V'/) - CHARACTER(1) iaxis_char,item_char - INTEGER iaxis,item,iion,ion_tmp - - read(datastring(6:),*)ion_tmp - do iion=1,nitot_max_forces - if(ion_tmp==iion)then -! Generate the number of atoms for which forces are calculated. - if(iion>nitot_forces)nitot_forces=iion - do iaxis=1,3 - iaxis_char=axis(iaxis) - if(datastring(5:5)==iaxis_char)then -! Generate the number of axis for which forces are calculated. - if(iaxis>naxis_forces)naxis_forces=iaxis - do item=1,22 - item_char=atem(item) - if(datastring(4:4)==item_char)then -! Generate the number of items. - if(item>nitem_forces)nitem_forces=item - call check_tag_free(tag_forces(item,iaxis,iion)) -! Generate label for forces - tag_forces(item,iaxis,iion)=i - endif ! item_char - enddo ! item - endif ! datastring=iaxis_char - enddo ! iaxis - endif ! iion - enddo ! iion - - END SUBROUTINE generate_tag_forces - - - SUBROUTINE construct_write_forces(startline,Nstudy,block_length) -!----------------------------------------------------------------! -! This routine calculates VMC/DMC forces from available data and ! -! performs a reblocking analysis. ! -! ! -! AB 11.2007 ! -!----------------------------------------------------------------! - IMPLICIT NONE - INTEGER,INTENT(in) :: startline,Nstudy,block_length - CHARACTER(20) plotname - INTEGER i,n,ialloc,nthird,nthirdstart,nthirdstop,ierr - DOUBLE PRECISION etot,std_err,delta_std_err,etot_dmc,etot_dmc_SE - DOUBLE PRECISION av,var,skew,kurt,max_val,min_val - LOGICAL forces_reblock,ltemp - -! Do we want to reblock forces? - forces_reblock=.false. - do - write(6,*)'Forces data are detected. When you like to reblock the forces da& - &ta with the' - write(6,*)'same block length as the total energy, choose F. When you like t& - &o investigate' - write(6,*)'reblocked forces error bars, choose T and use gnuplot to look at& - & error bars.' - write(6,*)'Choose F or T:' - read(5,*,iostat=ierr)forces_reblock - if(ierr/=0)forces_reblock=.false. - if((forces_reblock).or.(.not.forces_reblock))then - exit - else - write(6,*)'Please try again. Choose T or F.' - endif - enddo ! choice loop - - - write(6,*)'FORCES COMPONENTS WITH REBLOCKED ERROR BARS' - write(6,*)'===========================================' - -!------------------ reblock VMC forces ------------------------- - - if(trim(qmc_method)=='VMC')then - - allocate(forces_array(Nlines,11),stat=ialloc) - if(ialloc/=0)then - write(6,*)'Allocation problem (1).' - stop - endif - -! Need energy estimate - call reblock_unweighted(Nstudy,data_array(startline:Nlines,& - &tag_energy),1,etot,std_err,delta_std_err) - - inquire(file='DMC_energy',exist=ltemp) - if(ltemp)then - open(11,file='DMC_energy') - read(11,*)etot_dmc,etot_dmc_SE - close(11) - else - etot_dmc=0.d0 - endif - - do iion=1,nitot_forces - do iaxis=1,naxis_forces - write(6,*)'Forces on atom ',trim(i2s(iion)),' along axis ',& - &trim(i2s(iaxis)),' Mean (au) Err (au)' - -! Construct various VMC estimators for the forces - do n=1,Nlines -! 1. Total forces (d-loc) - forces_array(n,1)=data_array(n,tag_forces(9,iaxis,iion))& - &-2.d0*data_array(n,tag_forces(2,iaxis,iion))& - &+2.d0*etot*data_array(n,tag_forces(1,iaxis,iion)) -! HFT forces (d-loc) - forces_array(n,2)=data_array(n,tag_forces(9,iaxis,iion)) - if(data_array(n,tag_forces(10,iaxis,iion))/=0.d0)then -! HFT forces (p-loc) - forces_array(n,3)=data_array(n,tag_forces(10,iaxis,iion)) -! HFT forces (s-loc) - forces_array(n,4)=data_array(n,tag_forces(11,iaxis,iion)) - endif -! Wavefunction Pulay term - forces_array(n,5)=& - &-2.d0*data_array(n,tag_forces(2,iaxis,iion))& - &+2.d0*etot*data_array(n,tag_forces(1,iaxis,iion)) -! Pseudopotential Pulay term - forces_array(n,6)=-data_array(n,tag_forces(7,iaxis,iion))& - &+data_array(n,tag_forces(4,iaxis,iion)) -! 2. Total forces zero-variance corrected (class 1,d-loc) - forces_array(n,7)=data_array(n,tag_forces(9,iaxis,iion))& ! HFT - &-2.d0*data_array(n,tag_forces(2,iaxis,iion))& ! Pulay - &+2.d0*etot*data_array(n,tag_forces(1,iaxis,iion))& ! " - &-data_array(n,tag_forces(6,iaxis,iion))& ! -H Psi' - &-data_array(n,tag_forces(7,iaxis,iion))& ! " - &+data_array(n,tag_forces(3,iaxis,iion))& ! " - &+data_array(n,tag_forces(4,iaxis,iion)) ! " -! Zero-variance term - forces_array(n,8)=& - &-data_array(n,tag_forces(6,iaxis,iion))& !- H Psi' - &-data_array(n,tag_forces(7,iaxis,iion))& ! " - &+data_array(n,tag_forces(3,iaxis,iion))& ! " - &+data_array(n,tag_forces(4,iaxis,iion)) ! " - ! E_l Psi' term cancelled -! VMC nodal term, added to Total Forces (purHFT,purNT,d-loc) - if(etot_dmc/=0.d0)then - forces_array(n,9)=& ! - &-data_array(n,tag_forces(6,iaxis,iion))& ! -H Psi' - &-data_array(n,tag_forces(7,iaxis,iion))& ! " - &-data_array(n,tag_forces(2,iaxis,iion))& ! " - &+data_array(n,tag_forces(3,iaxis,iion))& ! " - &+data_array(n,tag_forces(4,iaxis,iion))& ! " - &+etot_dmc*data_array(n,tag_forces(1,iaxis,iion)) ! E*Psi' - else - forces_array(n,9)=0.d0 - endif - enddo ! Nlines - -! Write out the various estimators for the forces -! 1. Total forces - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,1),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'Total Force(dloc) :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,2),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'HFT Force(dloc) :',av,std_err - if(data_array(n,tag_forces(10,iaxis,iion))/=0.d0)then - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,3),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'HFT Force(ploc) :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,4),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'HFT Force(sloc) :',av,std_err - endif - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,5),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'Wavefunction Pulay term :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,6),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'Pseudopotential Pulay term :',av,std_err - -! 2. Total forces zero-variance corrected - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,7),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'Total Force+ZV(dloc) :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,8),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'Zero-variance term :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,9),& - &block_length,av,std_err,delta_std_err) - write(6,9) 'VMC NT(add to last DMC est) :',av,std_err - - if(forces_reblock)then - plotname='forces'//'.tot.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:& - &Nlines,1),plotname) - plotname='forces'//'.HFT.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:& - &Nlines,2),plotname) - plotname='forces'//'.totZV.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:& - &Nlines,7),plotname) - plotname='forces'//'.vmcNT.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:& - &Nlines,9),plotname) -! Analyse total forces by thirds if there is enough data. - if(Nstudy>=6)then - write(6,*) - write(6,*)'ANALYSIS OF TOTAL FORCES DATA BY THIRDS' - nthird=Nstudy/3 - write(6,*)' Data range Av forces Variance Maximum & - &Minimum' - do i=1,3 - nthirdstart=startline+(i-1)*nthird - nthirdstop=nthirdstart+nthird-1 - call compute_stats_unweighted(.false.,nthird, & - &forces_array(nthirdstart:nthirdstop,1),av,var,skew,kurt, & - &max_val,min_val) - write(6,'(" ",a16,4(" ",es12.5))')trim(i2s(nthirdstart)) & - &//'->'//trim(i2s(nthirdstop)),av,var,max_val,min_val - enddo ! i - else - write(6,*)'Not enough data to analyse by thirds: need at least 6 points.' - endif ! Enough data? - endif ! reblock - write(6,*) - - enddo ! iaxis - enddo ! iion - write(6,*) "The last estimator 'VMC NT' is zero (not used), unless a 'DMC_en& - &ergy' file is " - write(6,*) "provided in the working directory during the reblocking process & - &which contains" - write(6,*) "the DMC energy followed by its error bar. After this estimator V& - &MC NT is " - write(6,*) "calculated by the reblocking routine, please add it by hand to t& - &he 'Total" - write(6,*) "Force(purHFT/NT,dloc)' estimator to obtain another force estimat& - &or." - write(6,*) - endif ! VMC - -!------------------- reblock DMC forces --------------------- - - if(trim(qmc_method)=='DMC')then - allocate(forces_array(Nlines,10),stat=ialloc) - if(ialloc/=0)then - write(6,*)'Allocation problem (1).' - stop - endif - -! Need energy estimate - call reblock_weighted(Nstudy,data_array(startline:Nlines,tag_energy),& - &data_array(startline:Nlines,tag_weight),1,etot,std_err,delta_std_err) - -! Construct various DMC estimators for the forces - do iion=1,nitot_forces - do iaxis=1,naxis_forces - write(6,*)'Forces on atom ',trim(i2s(iion)),' along axis ',& - &trim(i2s(iaxis)),' Mean (au) Err (au)' - - do n=1,Nlines -! 1. Pure total forces (mixNT,dloc) - forces_array(n,1)=data_array(n,tag_forces(20,iaxis,iion))& ! HFT - &-data_array(n,tag_forces(18,iaxis,iion))& ! PPT - &+data_array(n,tag_forces(15,iaxis,iion))& ! " - &-2.d0*(data_array(n,tag_forces(6,iaxis,iion))& ! NT - & +data_array(n,tag_forces(8,iaxis,iion))& ! " - & +data_array(n,tag_forces(2,iaxis,iion))& ! " - & -data_array(n,tag_forces(3,iaxis,iion))& ! " - & -data_array(n,tag_forces(5,iaxis,iion))& ! " - & -etot*data_array(n,tag_forces(1,iaxis,iion))) ! " -! Pure total forces (purNT,dloc) - forces_array(n,2)=data_array(n,tag_forces(20,iaxis,iion))& ! HFT - &-data_array(n,tag_forces(18,iaxis,iion))& ! PPT - &+data_array(n,tag_forces(15,iaxis,iion))& ! " - &-(data_array(n,tag_forces(17,iaxis,iion))& ! NT - & +data_array(n,tag_forces(19,iaxis,iion))& ! " - & +data_array(n,tag_forces(13,iaxis,iion))& ! " - & -data_array(n,tag_forces(14,iaxis,iion))& ! " - & -data_array(n,tag_forces(16,iaxis,iion))& ! " - & -etot*data_array(n,tag_forces(12,iaxis,iion))) ! " -! Pure HFT forces (d-loc) - forces_array(n,3)=data_array(n,tag_forces(20,iaxis,iion)) - if(data_array(n,tag_forces(10,iaxis,iion))/=0.d0)then -! Pure HFT forces (p-loc) - forces_array(n,4)=data_array(n,tag_forces(21,iaxis,iion)) -! Pure HFT forces (s-loc) - forces_array(n,5)=data_array(n,tag_forces(22,iaxis,iion)) - endif -! Nodal term (mix) - forces_array(n,6)=& - &-2.d0*(data_array(n,tag_forces(6,iaxis,iion))& ! NT - & +data_array(n,tag_forces(8,iaxis,iion))& ! " - & +data_array(n,tag_forces(2,iaxis,iion))& ! " - & -data_array(n,tag_forces(3,iaxis,iion))& ! " - & -data_array(n,tag_forces(5,iaxis,iion))& ! " - & -etot*data_array(n,tag_forces(1,iaxis,iion))) ! " -! Nodal term (pur) - forces_array(n,7)=& - &-(data_array(n,tag_forces(17,iaxis,iion))& ! NT - & +data_array(n,tag_forces(19,iaxis,iion))& ! " - & +data_array(n,tag_forces(13,iaxis,iion))& ! " - & -data_array(n,tag_forces(14,iaxis,iion))& ! " - & -data_array(n,tag_forces(16,iaxis,iion))& ! " - & -etot*data_array(n,tag_forces(12,iaxis,iion))) ! " -! Pseudopotential Pulay term: Psi^(-1)WPsi'-Psi^(-1)WPsi Psi'/Psi - forces_array(n,8)=-data_array(n,tag_forces(18,iaxis,iion))& - &+data_array(n,tag_forces(15,iaxis,iion)) -! 2. Mixed total forces (d-loc) - forces_array(n,9)=data_array(n,tag_forces(9,iaxis,iion))& ! HFT - & -data_array(n,tag_forces(7,iaxis,iion))& ! PPT - & +data_array(n,tag_forces(4,iaxis,iion))& ! " - & -data_array(n,tag_forces(6,iaxis,iion))& ! NT - & -data_array(n,tag_forces(8,iaxis,iion))& ! " - & -2.d0*data_array(n,tag_forces(2,iaxis,iion))& ! " - & +data_array(n,tag_forces(3,iaxis,iion))& ! " - & +data_array(n,tag_forces(5,iaxis,iion))& ! " - & +2.d0*etot*data_array(n,tag_forces(1,iaxis,iion)) ! " - -! Mixed HFT forces (d-loc) - forces_array(n,10)=data_array(n,tag_forces(9,iaxis,iion)) - enddo ! Nlines - -! Write out forces - if(tag_weight>0.and.use_weights)then -! 1. Pure total forces - call reblock_weighted(Nstudy,forces_array(startline:Nlines,1),& - &data_array(startline:Nlines,tag_weight),& - &block_length,av,std_err,delta_std_err) - write(6,9)'Total Force(purHFT,mixNT,dloc) :',av,std_err - call reblock_weighted(Nstudy,forces_array(startline:Nlines,2),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'Total Force(purHFT,purNT,dloc) :',av,std_err - call reblock_weighted(Nstudy,forces_array(startline:Nlines,3),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'HFT Force(pur,dloc) :',av,std_err - if(data_array(n,tag_forces(10,iaxis,iion))/=0.d0)then - call reblock_weighted(Nstudy,forces_array(startline:Nlines,4),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'HFT Force(pur,ploc) :',av,std_err - call reblock_weighted(Nstudy,forces_array(startline:Nlines,5),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'HFT Force(pur,sloc) :',av,std_err - endif - call reblock_weighted(Nstudy,forces_array(startline:Nlines,6),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'Nodal Term(mix) :', av,std_err - call reblock_weighted(Nstudy,forces_array(startline:Nlines,7),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'Nodal Term(pur) :', av,std_err - call reblock_weighted(Nstudy,forces_array(startline:Nlines,8),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'Pseudopot. Pulay Term(pur) :',av,std_err -! Mixed total forces - call reblock_weighted(Nstudy,forces_array(startline:Nlines,9),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'Total Force(mix,dloc) :',av,std_err - call reblock_weighted(Nstudy,forces_array(startline:Nlines,10),& - &data_array(startline:Nlines,tag_weight),block_length,av,std_err,& - &delta_std_err) - write(6,9)'HFT Force(mix,dloc) :',av,std_err - write(6,*) - -! Reblock forces and write out into files - if(forces_reblock)then - plotname='forces'//'.totpur.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:Nlines,1),& - &plotname,data_array(startline:Nlines,tag_weight)) - plotname='forces'//'.HFTpur.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:Nlines,3),& - &plotname,data_array(startline:Nlines,tag_weight)) - plotname='forces'//'.totmix'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:Nlines,9),& - &plotname,data_array(startline:Nlines,tag_weight)) - plotname='forces'//'.HFT'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:Nlines,10),& - &plotname,data_array(startline:Nlines,tag_weight)) - endif ! forces_reblock - - else ! have weights -! 1. Pure total forces - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,1),& - &block_length,av,std_err,delta_std_err) - write(6,9)'Total Force(purHFT,mixHFT,dloc) :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,2),& - &block_length,av,std_err,delta_std_err) - write(6,9)'Total Force(purHFT,purHFT,dloc) :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,3),& - &block_length,av,std_err,delta_std_err) - write(6,9)'HFT Force(pur,dloc) :',av,std_err - if(data_array(n,tag_forces(10,iaxis,iion))/=0.d0)then - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,4),& - &block_length,av,std_err,delta_std_err) - write(6,9)'HFT Force(pur,ploc) :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,5),& - &block_length,av,std_err,delta_std_err) - write(6,9)'HFT Force(pur,sloc) :',av,std_err - endif - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,6),& - &block_length,av,std_err,delta_std_err) - write(6,9)'Nodal Term(mix) :', av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,7),& - &block_length,av,std_err,delta_std_err) - write(6,9)'Nodal Term(pur) :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,8),& - &block_length,av,std_err,delta_std_err) - write(6,9)'Pseudopot. Pulay Term(pur) :',av,std_err -! Mixed total forces - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,9),& - &block_length,av,std_err,delta_std_err) - write(6,9)'Total Force(mix,dloc) :',av,std_err - call reblock_unweighted(Nstudy,forces_array(startline:Nlines,10),& - &block_length,av,std_err,delta_std_err) - write(6,9)'HFT Force(mix,dloc) :',av,std_err - write(6,*) - -! Reblock forces and write out into files - if(forces_reblock)then - plotname='forces'//'.totpur.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:Nlines,1),& - &plotname) - plotname='forces'//'.HFTpur.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:Nlines,3),& - &plotname) - plotname='forces'//'.totmix.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:Nlines,9),& - &plotname) - plotname='forces'//'.HFTmix.'//trim(i2s(iaxis))//trim(i2s(iion))//& - &'.plot' - call reblock_forces_analysis(Nstudy,forces_array(startline:Nlines,10),& - &plotname) - endif - endif ! have weights - - enddo ! naxis - enddo ! nitot_forces - - endif ! DMC - -9 format(" ",a32,2(" ",f20.14)) - - END SUBROUTINE construct_write_forces - - -END MODULE analysis - - -PROGRAM analyse_qmc -!---------------------------! -! Main program starts here. ! -!---------------------------! - USE analysis, ONLY : filename,read_data,check_data,compute_stats - IMPLICIT NONE - LOGICAL vmc,dmc - - write(6,*) - write(6,*)'O---------O' - write(6,*)'| REBLOCK |' - write(6,*)'O---------O' - write(6,*) - -! What files are present? - inquire(file='vmc.hist',exist=vmc) - inquire(file='dmc.hist',exist=dmc) - - if(.not.(vmc.or.dmc))then - write(6,*)'Sorry, there are no vmc.hist or dmc.hist files to analyse.' - stop - endif ! No hist files found. - -! Sort out which file to analyse if more than one possibility exists. - if(dmc)then - filename='dmc.hist' - elseif(vmc)then - filename='vmc.hist' - else - write(6,*)'Bug.' - stop - endif - write(6,*)'Data in '//trim(filename)//' will be analysed.' - write(6,*) - -! Read in data from the file. - call read_data(dmc) - -! Check the data for inconsistencies and get units etc. - call check_data - -! Analyse the data. - call compute_stats - - write(6,*)'Program finished.' - write(6,*) - -END PROGRAM analyse_qmc