main.f90 Source File

#include "precompilerdefinitions"
program samples_from_md
!!{!src/samples_from_md/manual.md!}
use konstanter, only: flyt,lo_huge,lo_tol
use helpers, only: open_file,lo_mean,lo_stddev,lo_random_int,tochar,lo_seed_random_numbers,lo_progressbar_init,&
                   lo_progressbar,lo_shuffle_int_array
use type_crystalstructure, only: lo_crystalstructure
use type_mdsim, only: lo_mdsim
!use dump_data
use options, only: lo_opts

implicit none
type(lo_opts) :: opts
type(lo_mdsim) :: sim
type(lo_crystalstructure) :: p

integer, dimension(:), allocatable  :: sample, indices, testsample
integer :: i,j,k,ii,jj,mindist,u
real(flyt) :: f0,f1,f2,t0,t1
real(flyt) :: mean_ep,mean_ek,sigma_ep,sigma_ek
real(flyt) :: mp0,mp1,mk0,mk1,sp0,sp1,sk0,sk1
character(len=80) :: fname

! Get options
call opts%parse()
! Seed random numbers
call lo_seed_random_numbers()
! Get simulation
call sim%read_from_file()

! List of all the samples
allocate(indices(sim%nt))
do i=1,sim%nt
    indices(i)=i
enddo

! Get a first sample
allocate(sample(opts%n))
allocate(testsample(opts%n))
call random_number(f0)
sample(1)=floor(f0*sim%nt)+1

! Get the minimum distance between two samples
mindist=sim%nt/opts%n/5
write(*,*) 'minimum distance',mindist

! Choose the rest
do i=2,opts%n
    j=select_ok_random_sample(indices,sample(1:i-1),mindist)
    sample(i)=j 
enddo

! Values for the entire simulation
mean_ep=lo_mean(sim%stat%ep)
mean_ek=lo_mean(sim%stat%ek)
sigma_ep=lo_stddev(sim%stat%ep)
sigma_ek=lo_stddev(sim%stat%ek)

! Starting values for monte-carlo run
mp0=lo_mean(sim%stat%ep(sample))
mk0=lo_mean(sim%stat%ek(sample))
sp0=lo_stddev(sim%stat%ep(sample))
sk0=lo_stddev(sim%stat%ek(sample))
mp0=abs(mp0-mean_ep)
mk0=abs(mk0-mean_ek)
sp0=abs(sp0-sigma_ep)
sk0=abs(sk0-sigma_ek)

t0=mp0+mk0+sp0+sk0

write(*,*) 'Starting Monte-Carlo, distances:'
write(*,"(9X,4(1X,A18))") '|<Ep>-<Ep0>|','|<Ek>-<Ek0>|','|S(Ep)-S(Ep0)|','|S(Ek)-S(Ek0)|'
write(*,"(1X,I8,4(1X,F18.6))") 0,mp0,mk0,sp0,sk0

do i=1,opts%maxiter
    ! try a new configuration
    testsample=sample
    ! CHANGE PLACES
    k=lo_random_int(opts%n)
    j=select_ok_random_sample(indices,sample,mindist)
    testsample(k)=j
    
    ! Measure distances
    mp1=lo_mean(sim%stat%ep(testsample))
    mk1=lo_mean(sim%stat%ek(testsample))
    sp1=lo_stddev(sim%stat%ep(testsample))
    sk1=lo_stddev(sim%stat%ek(testsample))
    mp1=abs(mp1-mean_ep)
    mk1=abs(mk1-mean_ek)
    sp1=abs(sp1-sigma_ep)
    sk1=abs(sk1-sigma_ek)

    t1=mp1+mk1+sp1+sk1
    ! MC check
    call random_number(f2)
    if ( exp(-(t1-t0)/opts%temp) .gt. f2 ) then
        sample=testsample
        mp0=mp1
        mk0=mk1
        sp0=sp1
        sk0=sk1
        t0=t1
        write(*,"(1X,I8,4(1X,F18.6))") i,mp0,mk0,sp0,sk0
        if ( t0 .lt. lo_tol ) exit
    endif
enddo

write(*,*) ' '
write(*,*) ' Results of Monte-Carlo run '
write(*,"('                 ',4(1X,A18))") '<Ep>','<Ek>','Sigma(Ep)','Sigma(Ek)'
write(*,"('Full simulation:',4(1X,F18.6))") &
        lo_mean(sim%stat%ep),lo_mean(sim%stat%ek),lo_stddev(sim%stat%ep),lo_stddev(sim%stat%ek)
write(*,"('        Samples:',4(1X,F18.6))") &
        lo_mean(sim%stat%ep(sample)),lo_mean(sim%stat%ek(sample)),lo_stddev(sim%stat%ep(sample)),lo_stddev(sim%stat%ek(sample))

! Find the minimum, maximum and mean distance between two samples
ii=0
jj=1000000000
f0=0.0_flyt
do i=1,opts%n
    f1=lo_huge
    do j=1,opts%n
        if ( i .ne. j ) then
            jj=min(abs((sample(i)-sample(j))),jj)
            f1=min(abs((sample(i)-sample(j))*1.0_flyt),f1)
        endif
    enddo
    ii=max(floor(f1),ii)
    f0=f0+f1
enddo
f0=f0/opts%n
write(*,*) '  Max distance:',ii
write(*,*) '  Min distance:',jj
write(*,*) ' Mean distance:',int(f0)

! write the samples
p=sim%crystalstructure
call lo_shuffle_int_array(sample)
do i=1,opts%n
    p%r=sim%r(:,:,sample(i))
    p%v=sim%v(:,:,sample(i))
    p%info%title='sample '//tochar(sample(i),7)
    fname='sample'//tochar(i,5)
    select case(opts%output_format)
    case(1) ! VASP
        call p%writetofile(fname,opts%output_format,.true.)
    case(2) ! Abinit                                       
        call p%writetofile(fname,opts%output_format,.true.)
    case default
        write(*,*) 'Output format not fixe yet'
        stop
    end select
enddo

! and the energies and stuff of the samples
u=open_file('out','outfile.stat_sample')
    do j=1,opts%n
        i=sample(j)
        write(u,"(I6,' ',13(2X,EN16.7))") &
        i,i*sim%ts,sim%stat%et(i),sim%stat%ep(i),sim%stat%ek(i),sim%stat%t(i),sim%stat%p(i),sim%stat%stress(:,i)
    enddo
close(u)

contains

!> choose a sample, but not too stupidly
integer function select_ok_random_sample(indices,forbidden,mindist)
    integer, dimension(:), intent(in) :: indices, forbidden
    integer, intent(in) :: mindist
    !
    integer :: i,ii,j,l,nt,nf
    
    nf=size(forbidden,1)
    nt=size(indices,1)
    do ii=1,100
        i=lo_random_int(nt)
        l=0
        do j=1,nf
            if ( abs(indices(i)-indices(forbidden(j))) .le. mindist ) l=l+1
        enddo
        if ( l .eq. 0 ) exit
    enddo
    select_ok_random_sample=i
end function

end program