mirror of
https://github.com/pfloos/quack
synced 2024-08-25 21:51:46 +02:00
2765 lines
100 KiB
Fortran
2765 lines
100 KiB
Fortran
|
!--------------------------------------------------------------------!
|
||
|
! 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)))
|
||
|
|
||
|
! <O>, <O>**2, <O**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<tol)return
|
||
|
invOvar=1.d0/Ovar
|
||
|
|
||
|
! Autocorrelation for i<=cut-off -> 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)econst_hist=0.d0
|
||
|
if(econst_is_const)then
|
||
|
if(abs(constant_energy-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<Nequil)then
|
||
|
write(6,*)'Warning: equilibration data will be included in the &
|
||
|
&statistics analysis.'
|
||
|
write(6,*)
|
||
|
endif
|
||
|
else
|
||
|
! Don't skip any lines for VMC.
|
||
|
Nskip=0
|
||
|
endif ! QMC method
|
||
|
write(6,*)trim(i2s(Nlines-Nskip))//' lines of data will be analysed.'
|
||
|
write(6,*)'These data start on line '//trim(i2s(Nskip+1)) &
|
||
|
&//' and end on line '//trim(i2s(Nlines))//'.'
|
||
|
write(6,*)
|
||
|
|
||
|
! Find out how input data is scaled w.r.t. the simulation cell
|
||
|
if(trim(atom_basis_type)=='none')then
|
||
|
nparts_per_simcell=netot
|
||
|
elseif(isperiodic)then
|
||
|
nparts_per_simcell=npcells
|
||
|
else
|
||
|
nparts_per_simcell=1
|
||
|
endif
|
||
|
|
||
|
! Find out what units are to be used.
|
||
|
if(trim(atom_basis_type)=='none')then
|
||
|
! Electron gas.
|
||
|
do
|
||
|
write(6,*)'Your data are for an electron(-hole) system.'
|
||
|
write(6,*)'Please select units for your energy data.'
|
||
|
write(6,*)'Choose one of: (1) au per particle; &
|
||
|
&(2) eV per particle.'
|
||
|
read(5,*,iostat=ierr)units_choice
|
||
|
if(ierr/=0)units_choice=-1
|
||
|
if(units_choice>=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
|