!8/22/01 ! lame6.hpf ! STATUS: incomplete, workable ! -must find cause of "segmentation fault" ! -re-modifying the correlation to be CORRECTLY periodic- corrbf2 ! and corrln have to be completed ! QUESTIONS: ! -is the periodic correlation function correct in the program ! -related: is it okay to take only the log of absolute values in corrln ! Salman says we need to normalize the correlation function. Why isn't ! this being done already? Although it is clear that the corrln will never ! operate correctly if we're getting negative values for corr., since i take ! the abs value before taking the log !------------------------------------------------------------------------------ ! lame6.hpf ! Matt Holmes, LANL T-8 ! Based on the program qes.hpf written by Salman Habib. ! This program evolves a field on a lattice using the lame potential. ! The lame potential is periodic, formed from the square of an elliptical ! function. Its special properties of interest to us is that the equivalent ! quantum mechanical schrodinger equation (obtained via the transfer integral ! method) is exactly solvable. This will allow us to compare the numerical ! correlation length and pdf calculated by this program with the analytic ! results. ! Log of modifications from the previous qes7.hpf: ! A) Include subroutine to calculate sn(x,k),cn(x,k),dn(x,k). Subroutine ! comes from Thompson. ! B) Modify subroutine "pot" to print out the lame potential (uses sn). ! C) Modify subroutine "rk2" to implement the lame force (uses sn,dn,cn). ! D) Modify subroutine "pot" to print out force also (if desired). ! E) Calculate periodic correlations and correlation length ! F) Take log only of absolute values in corrln ! G) Subroutine perpdf will periodify the pdfsym, output pdf_sum ! H) Calculate periodic correlation using a cosine correlator ! I) Increase the number of output files in order to "restart" the simulation. ! ----------------------------------------------------------------------------- ! LIST OF OUTPUT FILES ("!" signifies that file output is commented out): ! NOTE: We care most about 85 and 96, with the pdf and the correlation length. ! fort.24 x,sn,cn,dn ! fort.25 x,potential ! fort.26 x,force ! fort.60 itime, penergy, tenergy, kenergy ! fort.64 EQUIL: xlatticepoint, averaged field-field correlation, ! fld^2-fld^2 correlation ! fort.65 FINAL: xlatticepoint, averaged field-field correlation, ! fld^2-fld^2 correlation ! fort.80 EQUIL: latticebin, pdf, pdfsym ! fort.81 FINAL: latticebin, pdf, pdfsym ! fort.85 FINAL: phi,perpdf (for "periodic pdf") ! fort.90 grnd state energy values ! fort.96 xlatticepoint, scaled log field-field correlation, same for squared corr. ! fort.195 every data set of the field-field correlations, at each checkpoint ! fort.98 x,v ! fort.99 x,phi ! fort.311 phi ! fort.312 v ! fort.313 sumc1 ! fort.314 sumc2 ! fort.315 prevtime ! fort.316 pdf PROGRAM lame !-----Langevin solver for double well QES field theory. !-----Salman Habib, April 27, 1997. (uses Grant's stochastic RK2) implicit none !-----DEFINE INPUT VARIABLES------------------------------------------------- integer :: nx,np,npdf,pdstart,integer,pdfint integer :: interval,tottime,newdata,stime,naves real :: eps,theta,etabar,length,xxl,n,k print *, "program is running..." print *, "opening input file now" open(33,file='lame6.in') read(33,*)newdata !A logial saying whether to continue an old sim. read(33,*)theta !Temperature read(33,*)etabar !Viscosity coefficient read(33,*)nx !Number of spatial lattice points read(33,*)np !Number of correlation function lattice points. read(33,*)npdf !Number of histrogramming bins for the dist function read(33,*)eps !Delta t read(33,*)length !Delta x read(33,*)pdfint !Normalizing factor for the dist function read(33,*)pdstart !Time at which to begin taking data (after equil) read(33,*)interval!Amount of time to evolve between diagnostic checks read(33,*)stime !Start time read(33,*)tottime !Total simulation time read(33,*)xxl !Used to determine range of potential to plot read(33,*)n !A parameter in the lame potential read(33,*)k !The "modulus" for the elliptical function in the ! lame potential read(33,*)naves !The number of images of the shifted pdf to average close(33) print *, "Read input file." !-----CALL SOLVE (Master Program)-------------------------------------------- call solve(newdata,theta,etabar,nx,np,npdf,eps,length, & pdfint,pdstart,interval,stime,tottime,n,k,xxl,naves) end PROGRAM lame !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !-----SUBROUTINE SOLVE (Master Program)------------------------------------- SUBROUTINE solve(newdata,theta,etabar,nx,np,npdf,eps,length, & pdfint,pdstart,interval,stime,tottime,n,k,xxl,naves) implicit none !-----FRONT END VARIABLES---------------------------------------------------- ! Dummy variables passed to subroutine: integer,intent(in) :: newdata,nx,np,npdf,pdfint,pdstart, & interval,stime,tottime,naves real,intent(in) :: theta,etabar,eps,length,n,k,xxl ! Local variables: integer :: i,itime,otime,ftime,prevtime real :: int,spacing integer,dimension(npdf) :: pdf,add !hpf$ distribute pdf(block) !hpf$ align (:) with pdf(:) :: add real,dimension(npdf) :: pdfsym, pdf_sum !hpf$ distribute pdfsym(block) !hpf$ align (:) with pdfsym(:) :: pdf_sum real,dimension(nx) :: phi,noise,v !hpf$ distribute phi(block) !hpf$ align (:) with phi(:) :: noise,v real,dimension(np) :: sumc1,sumc2 real,dimension(np) :: c1add,c2add,sumc1w,sumc2w !hpf$ distribute sumc1(block) !hpf$ align (:) with sumc1(:) :: sumc2,c1add,c2add,sumc1w,sumc2w !-------F90 RANDOM NUMBER STUFF----------------------------------------------- ! why the hell is this needed? we don't change the value of the seed, and ! we don't use this array anywhere else ! integer,dimension(:),allocatable :: seed integer seed_size call random_seed(SIZE=seed_size) print *, "seed_size is ", seed_size ! allocate (seed(seed_size)) ! call random_seed(GET=seed(1:seed_size)) ! if we want to set the seed, we should read "seed" from the input file, and ! use this next line ! random_seed(PUT=seed) ! where seed should be an integer larger than seed_size print *, "Allocated main arrays, set random seed size" !-------PRINT OUT POTENTIAL IF NEEDED----------------------------------------- print *, "Printing the potential to file 'fort.25'" call pot(n,k,xxl,spacing) !-------READ IN INPUT FILES IF NEEDED----------------------------------------- ! The information we need to continue the simulation is phi and v. ! The information we need to continue the diagnostics are sumc1,sumc2,pdf; ! also, we need to keep track of the number of data sets to average ! over; I stored this as "prevtime". if (newdata.eq.1) then print *, "Continuing old simulation, reading files." open(311,file='fort.311',status='unknown') read(311,*)phi close(311) open(312,file='fort.312',status='unknown') read(312,*)v close(312) open(313,file='fort.313',status='unknown') read(313,*)sumc1 close(313) open(314,file='fort.314',status='unknown') read(314,*)sumc2 close(314) open(315,file='fort.315',status='unknown') read(315,*)prevtime close(315) open(316,file='fort.316',status='unknown') read(316,*)pdf close(316) else if (newdata.eq.0) then !-------INITIALIZE PHI IF FILES ARE NOT READ IN------------------------------- phi=0 v=0 sumc1=0.0 sumc2=0.0 pdf=0.0 prevtime=0 end if !-------INITIALIZATIONS------------------------------------------------------- pdfsym=0 pdf_sum=0 otime=0 c1add=0 c2add=0 !-------PERFORM EQUILIBRATION ONLY FOR NEW SIMULATION if (newdata.eq.0) then print *, "Starting equilibration loop now" !-------LOOPING STARTS (INITIAL LOOP)----------------------------------------- do itime=stime,pdstart !-------2nd order RK Langevin step call random_number(noise) where(noise.gt.5.0/6.0)noise=sqrt(6.0*etabar*theta/length/eps) where(noise.le.1.0/6.0)noise=-sqrt(6.0*etabar*theta/length/eps) where(abs(noise).le.1.0)noise=0.0 call rk2(phi,v,noise,length,etabar,eps,n,k,nx) enddo print *, "Finished equil. loop, performing equil diagnostics." !-------Call diagnostics/analysis routines itime=pdstart print *, "Calling unpdf:" call unpdf(phi,add,nx,pdfint,npdf) pdf=pdf+add !UNCOMMENT for hpf: forall(i=1:npdf) pdfsym(i)=(pdf(i)+pdf(npdf+1-i))/2.0 !UNCOMMENT for f90: ! do i=1,npdf ! pdfsym(i)=(pdf(i)+pdf(npdf+1-i))/2.0 ! end do pdfsym=pdfsym/(1.0*(prevtime+1)) int=(sum(pdfsym))/(1.0*pdfint) pdfsym=pdfsym/int ! call energy(phi,v,length,n,k,nx,itime) ! call grndenergy(pdfsym,npdf,pdfint,itime,n,k,theta) print *, "Post equilibration, calling corrbf1" call corrbf1(phi,c1add,nx,np,spacing) sumc1=sumc1+c1add sumc1w=sumc1/(1.0*(prevtime+1)) !Average the correlations. print *, "Post equilibration, calling corrbf2" call corrbf2(phi,c2add,nx,np,spacing) sumc2=sumc2+c2add sumc2w=sumc2/(1.0*(prevtime+1)) !Average the correlations. print *, "Calling perpdf" call perpdf(npdf,pdfsym,pdf_sum,spacing,pdfint,naves) print *, "Finished equilibrium loop, performing diagnostics..." call printi(pdf,pdfsym,npdf,pdfint,np,sumc1w,sumc2w) end if print *, "starting main loop" !-------LOOPING STARTS (NESTED LOOP)------------------------------------------ !-------Outer Loop: Carry out diagnostics/analysis at end of each step do otime=1,(tottime-pdstart)/interval print *, "One iteration of main loop going now..." !-------Second Loop: Time step through without interruptions do ftime=1,interval call random_number(noise) where(noise.gt.5.0/6.0)noise=sqrt(6.0*etabar*theta/length/eps) where(noise.le.1.0/6.0)noise=-sqrt(6.0*etabar*theta/length/eps) where(abs(noise).le.1.0)noise=0.0 call rk2(phi,v,noise,length,etabar,eps,n,k,nx) end do !-------Call diagnostics/analysis routines itime=pdstart+otime*interval print *, "In main loop, calling unpdf" call unpdf(phi,add,nx,pdfint,npdf) pdf=pdf+add ! UNCOMMENT for hpf: forall(i=1:npdf) pdfsym(i)=(pdf(i)+pdf(npdf+1-i))/2.0 ! UNCOMMENT for f90: ! do i=1,npdf ! pdfsym(i)=(pdf(i)+pdf(npdf+1-i))/2.0 ! end do pdfsym=pdfsym/(1.0*(otime+prevtime+1)) !We have done otime diag. in !main loop, 1 diag. in equil. loop, and prevtime !diag. in prev simulation= otime+prevtime+1 int=(sum(pdfsym))/(1.0*pdfint) pdfsym=pdfsym/int ! call energy(phi,v,length,n,k,nx,itime) ! call grndenergy(pdfsym,npdf,pdfint,itime,n,k,theta) print *, "in main loop, calling corrbf1" call corrbf1(phi,c1add,nx,np,spacing) sumc1=sumc1+c1add print *, "in main loop, calling corrbf2" call corrbf2(phi,c2add,nx,np,spacing) sumc2=sumc2+c2add end do print *, "Finished main loop, performing final diagnostics..." !-------Final computations/outputs sumc1w=sumc1/(1.0*(otime+prevtime+1)) !Average the correlations. sumc2w=sumc2/(1.0*(otime+prevtime+1)) prevtime=prevtime+otime !Record the number of data sets we have !averaged over,to be used in continuing !the simulation. print *, "otime= ",otime print *, "prevtime= ",prevtime print *, "Calling perpdf..." call perpdf(npdf,pdfsym,pdf_sum,spacing,pdfint,naves) print *, "Calling corrln (fort.96)..." call corrln(sumc1w,sumc2w,length,np) print *, "Calling printj (fort.65, fort.81,fort.85)..." call printj(pdf,pdfsym,pdf_sum,npdf,pdfint,np,sumc1w,sumc2w) print *, "Calling printfinal (fort.98,99,311,312,313,314,315,316)..." call printfinal(phi,v,nx,np,npdf,sumc1,sumc2,pdf,prevtime) end SUBROUTINE solve !-------SUBROUTINE PRINTI----------------------------------------------------- SUBROUTINE printi(pdf,pdfsym,npdf,pdfint,np,sumc1w, & sumc2w) implicit none ! Dummy variables: integer, intent(in) :: pdfint,npdf,np integer, dimension(npdf) :: pdf !hpf$ distribute pdf(block) real, dimension(npdf) :: pdfsym,xxx !hpf$ distribute pdfsym(block) !hpf align (:) with pdfsym(:) :: xxx real, dimension(np) :: sumc1w,sumc2w !hpf$ distribute sumc1w(block) !hpf align (:) with sumc1w(:) :: sumc2w ! Local variables: integer :: i, nxm real :: xxmin xxmin=-0.50*(1.0*npdf-1.0)/(1.0*pdfint) !UNCOMMENT for hpf: forall(i=1:npdf) xxx(i)=xxmin + (i-1)/(1.0*pdfint) !UNCOMMENT for f90: ! do i=1,npdf ! xxx(i)=xxmin + (i-1)/(1.0*pdfint) ! end do do i=1,np nxm=i-1 write(64,*)nxm,sumc1w(i),sumc2w(i) enddo do i=1,npdf write(80,*)xxx(i),pdf(i),pdfsym(i) enddo end SUBROUTINE printi !-------SUBROUTINE PRINTJ----------------------------------------------------- SUBROUTINE printj(pdf,pdfsym,pdf_sum,npdf,pdfint,np,sumc1w, & sumc2w) implicit none ! Dummy variables: integer, intent(in) :: pdfint,npdf,np integer, dimension(npdf) :: pdf !hpf$ distribute pdf(block) real, dimension(npdf) :: pdfsym,pdf_sum,xxx !hpf$ distribute pdfsym(block) !hpf$ align (:) with pdfsym(:) :: xxx,pdf_sum real, dimension(np) :: sumc1w,sumc2w !hpf$ distribute sumc1w(block) !hpf$ align (:) with sumc1w(:) :: sumc2w ! Local variables: integer :: i, nxm real :: xxmin xxmin=-0.50*(1.0*npdf-1.0)/(1.0*pdfint) !UNCOMMENT for hpf: forall(i=1:npdf) xxx(i)=xxmin + (i-1)/(1.0*pdfint) !UNCOMMENT for f90: ! do i=1,npdf ! xxx(i)=xxmin + (i-1)/(1.0*pdfint) ! end do do i=1,np nxm=i-1 write(65,*)nxm,sumc1w(i),sumc2w(i) enddo do i=1,npdf write(81,*)xxx(i),pdf(i),pdfsym(i) enddo do i=1,npdf write(85,*)xxx(i),pdf_sum(i) enddo end SUBROUTINE printj !-------SUBROUTINE PRINTFINAL-------------------------------------------------- SUBROUTINE printfinal(phi,v,nx,np,npdf,sumc1,sumc2,pdf,prevtime) implicit none ! Dummy variables: integer, intent(in):: nx,np,npdf real, dimension(nx), intent(in) :: phi,v real,dimension(np),intent(in) :: sumc1,sumc2 integer,dimension(npdf),intent(in) :: pdf integer,intent(in) :: prevtime ! Local variables: integer :: i !hpf$ distribute phi(block) !hpf$ align (:) with phi(:) :: v !hpf$ distribute sumc1(block) !hpf$ align (:) with sumc1(:) :: sumc2 !hpf$ distribute pdf(block) ! This format is for plotting: do i=1,nx write(99,*)i,phi(i) write(98,*)i,v(i) enddo ! The following output files ("300"series) are for continuing another sim. write(311,*)phi write(312,*)v write(313,*)sumc1 write(314,*)sumc2 write(315,*)prevtime write(316,*)pdf end SUBROUTINE printfinal !-------SUBROUTINE GRNDENERGY------------------------------------------------- SUBROUTINE grndenergy(pdfsym,npdf,pdfint,itime,n,k,theta) implicit none ! Dummy variables: integer,intent(in) :: pdfint,npdf,itime real,intent(in) :: n,k,theta real,dimension(npdf),intent(in) :: pdfsym ! Local variables: integer :: i real :: int,xxmin,kep,vep,tep real,dimension(npdf) :: xxx,psip real,dimension(npdf) :: sn,cn,dn !Holds values for elliptic functions !hpf$ distribute pdfsym(block) !hpf$ align (:) with pdfsym(:) :: xxx,psip,sn,cn,dn xxmin=-0.50*(1.0*npdf-1.0)/(1.0*pdfint) !UNCOMMENT for hpf: forall(i=1:npdf) xxx(i)=xxmin + (i-1)/(1.0*pdfint) !UNCOMMENT for f90: ! do i=1,npdf ! xxx(i)=xxmin + (i-1)/(1.0*pdfint) ! end do psip=sqrt(pdfsym) kep=-sum(psip*(eoshift(psip,-1,0.0,1)+eoshift(psip,1,0.0,1)-2.0* & psip)) *theta*theta*pdfint*0.50 call sncndn(xxx,k,npdf,sn,cn,dn) vep=(sum(pdfsym*(n*(n+1)*k*k*sn*sn)))/(1.0*pdfint) tep=kep+vep write(90,*)itime,tep,vep,kep end SUBROUTINE grndenergy !-------SUBROUTINE CORRLN------------------------------------------------- SUBROUTINE corrln(sumc1w,sumc2w,length,np) implicit none ! Dummy variables: integer,intent(in) :: np real, intent(in) :: length real,dimension(np),intent(in) :: sumc1w,sumc2w ! Local variables: integer :: i real,dimension(np) :: logc1,logc2 real,dimension(np-1) :: gamma1,gamma2 !hpf$ distribute sumc1w(block) !hpf$ align (:) with sumc1w(:) :: sumc2w,logc1,logc2 !hpf$ distribute gamma1(block) !hpf$ align (:) with gamma1(:) :: gamma2 logc1=log(abs(sumc1w)) logc2=log(abs(sumc2w)) !UNCOMMENT for hpf: forall(i=1:np-1) gamma1(i)=(logc1(i+1)-logc1(i))/length !UNCOMMENT for hpf: forall(i=1:np-1) gamma2(i)=(logc2(i+1)-logc2(i))/length !UNCOMMENT for f90: ! do i=1,np-1 ! gamma1(i)=(logc1(i+1)-logc1(i))/length ! end do !UNCOMMENT for f90: ! do i=1,np-1 ! gamma2(i)=(logc2(i+1)-logc2(i))/length ! end do do i=1,np-1 write(96,*)i,gamma1(i),gamma2(i) enddo end SUBROUTINE corrln !-------SUBROUTINE POT-------------------------------------------------------- ! Our potential is the lame potential. It is periodic and defined with the ! elliptical function sn. The parameters for this potential are n and k. ! The force is the derivative of the potential, and uses the functions sn, dn, ! and cn. ! This function also calculates the period of the potential for use in the ! per_pdf subroutine. SUBROUTINE pot(n,k,xxl,spacing) implicit none ! Dummy variables: real,intent(in) :: n,k,xxl !xxl gives the range of the pot to print real, intent(inout) :: spacing !spacing is the period of the potential ! Local variables: integer :: i,count real :: current,value real :: dx,xmin real,dimension(1024) :: potential,xx real,dimension(1024) :: sn,cn,dn real,dimension(1024) :: force !hpf$ distribute potential(block) !hpf$ align (:) with potential(:) :: xx,sn,cn,dn !hpf$ align (:) with potential(:) :: force xmin=-0.50*xxl dx=xxl/(1.0*1024) !UNCOMMENT for hpf: forall(i=1:1024) xx(i)=xmin+(i-1)*dx call sncndn(xx,k,1024,sn,cn,dn) !UNCOMMENT for hpf: potential=n*(n+1)*k*k*sn*sn force=-2.0*n*(n+1)*k*k*sn*cn*dn !UNCOMMENT for f90: ! do i=1,1024 ! xx(i)=xmin+(i-1)*dx ! end do !UNCOMMENT for f90: ! do i=1,1024 ! potential(i)=n*(n+1)*k*k*sn(i)*sn(i) ! force(i)=-2.0*n*(n+1)*k*k*sn(i)*cn(i)*dn(i) ! end do ! Here we calculate the spacing between minima of the potential, to be used ! to periodify the pdf. We do it with the force, by starting at 0 and finding ! the next zero of the force. Then we double it to obtain the period of the pot. count=512 !Start at the middle of the array, where xx(512)=0 current=0 value=0 spacing=0 do count=count+1 !will be used as an array indice (start at second !location in the array) current=force(count) !Get next value in the array if(current > 0) exit end do !checks: print *, "count is", count print *, "the first value greater than zero is", current print *, "the just previous value was", force(count-1) print *, "we are averaging", xx(count-1), " and", xx(count) value=(xx(count-1)+xx(count))/2 !use the average of the location just !before the zero and just after. spacing=2*value !The location of the first zero crossing is only !half the period of the potential. print *, "value is", value print *, "spacing is", spacing do i=1,1024 write(25,*)xx(i),potential(i) write(26,*)xx(i),force(i) write(24,*)xx(i),sn(i),cn(i),dn(i) !Write out the elliptic functions enddo end SUBROUTINE pot !-------SUBROUTINE CORRBF1---------------------------------------------------- SUBROUTINE corrbf1(phi,c1add,nx,np,spacing) ! Calculates a periodic field-field correlation function: ! C(u)=<|cos((pi/2)(1/min_location)(u(x)-u(0)))|> ! For this to be appropriately periodic, the values of the field must be ! scaled by pi/2; this is being achieved within the correlation calculation ! itself. implicit none ! Dummy variables: integer,intent(in) :: nx,np real,intent(in) :: spacing real,dimension(np),intent(inout) :: c1add real,dimension(nx),intent(in) :: phi !hpf$ distribute phi(block) ! Local variables: real,parameter :: pi=3.1415926, piby2=pi/2.0 integer :: nxm,i !nxm is used as a dummy index when writing to the file, see below real :: phiavg,phi2avg real,dimension(np) :: cphi !hpf$ distribute c1add(block) !hpf$ align (:) with c1add(:) :: cphi cphi=0 phiavg=sum(phi)/(1.0*nx) !averages will be used for normalization phi2avg=sum(phi*phi)/(1.0*nx) do i=1,np ! Calculate the correlation cphi(i)=sum(abs(cos((piby2/spacing)*(cshift(phi,i-1)-phi)))) $ /(1.0*nx) enddo ! NOT SURE WHAT NORMALIZATION TO USE NOW? ! c1add=(cphi-phiavg*phiavg)/(phi2avg-phiavg*phiavg) !a normalization c1add=cphi do i=1,np nxm=i-1 write(195,*)nxm,c1add(i) enddo end SUBROUTINE corrbf1 !-------SUBROUTINE CORRBF2---------------------------------------------------- SUBROUTINE corrbf2(phi,c2add,nx,np,spacing) implicit none ! Dummy variables: integer,intent(in) :: nx,np real,intent(in) :: spacing real,dimension(np),intent(inout) :: c2add real,dimension(nx),intent(in) :: phi !hpf$ distribute phi(block) ! Local variables: real,parameter :: pi=3.1415926, piby2=pi/2.0 integer :: nxm,i real :: phi2avg real,dimension(np) :: cphi2r,cphi2 !hpf$ distribute c2add(block) !hpf$ align (:) with c2add(:) :: cphi2,cphi2r cphi2=0.0 phi2avg=sum(phi*phi)/(1.0*nx) ! INCOMPLETE: how do we alter this to make it appropriately periodic? do i=1,np cphi2(i)=sum(abs(cos((piby2/spacing)*(cshift(phi,i-1)* $ cshift(phi,i-1)-phi*phi))))/(1.0*nx) enddo !A NORMALIZATION I DON'T UNDERSTAND: ! cphi2r=cphi2-phi2avg**2 c2add=cphi2r-phi2avg*phi2avg ! do i=1,np ! nxm=i-1 ! write(196,*)nxm,c2add(i) ! enddo end SUBROUTINE corrbf2 !-------SUBROUTINE UNPDF------------------------------------------------------ SUBROUTINE unpdf(phi,add,nx,pdfint,npdf) implicit none ! Dummy variables: integer,intent(in) :: nx,pdfint,npdf integer,dimension(npdf),intent(inout) :: add !hpf$ distribute add(block) real,dimension(nx) :: phi ! Local variables: integer :: i,nco logical,dimension(nx) :: mmsskk !hpf$ distribute mmsskk(block) real,dimension(nx) :: iphi !hpf$ distribute phi(block) !hpf$ align (:) with phi(:) :: iphi iphi=1.0*pdfint*phi do i=-(npdf-1)/2,(npdf-1)/2 mmsskk=.false. where(i.le.iphi.and.iphi.lt.i+1)mmsskk=.true. nco=count(mmsskk) add((npdf-1)/2+1+i)=nco enddo end SUBROUTINE unpdf !-------SUBROUTINE ENERGY----------------------------------------------------- SUBROUTINE energy(phi,v,length,n,k,nx,itime) implicit none ! Dummy variables: integer,intent(in) :: nx,itime real,intent(in) :: length,n,k real,dimension(nx),intent(in) :: phi !hpf$ distribute phi(block) ! Local variables: real :: kenergy,penergy,tenergy real,dimension(nx) :: gr,v real,dimension(nx) :: sn,cn,dn !hpf$ align (:) with phi(:) :: gr,v,sn,cn,dn gr=(cshift(phi,1,1)-phi)/length call sncndn(phi,k,nx,sn,cn,dn) penergy=length*sum(.50*gr*gr+n*(n+1)*k*k*sn*sn) kenergy=.50*length*sum(v*v) tenergy=kenergy+penergy penergy=penergy/(nx*length) kenergy=kenergy/(nx*length) tenergy=tenergy/(nx*length) write(60,*)itime,tenergy,penergy,kenergy end SUBROUTINE energy !-------SUBROUTINE RK2-------------------------------------------------------- SUBROUTINE rk2(phi,v,noise,length,etabar,eps,n,k,nx) implicit none ! Dummy variables: integer,intent(in) :: nx real,intent(in) :: length,etabar,eps,n,k real,dimension(nx),intent(in) :: noise real,dimension(nx),intent(inout)::phi,v ! Local variables: real,dimension(nx) :: sn,cn,dn !arrays to hold the values of elliptic !functions for each value of the phi array real,dimension(nx) ::pibefore,lapphi, phibefore,tv1,tv2,tphi1 !hpf$ align (:) with phi(:) :: sn,cn,dn !hpf$ align (:) with phi(:) :: pibefore,phibefore,lapphi,v,noise,tv1,tv2,tphi1 pibefore=v phibefore=phi lapphi=(cshift(phi,-1)+cshift(phi,1)-2.0*phi)/length/length ! Derivative of the potential is: ! dV/du = 2 n(n+1)k^2 cn(u,k)dn(u,k)sn(u,k) call sncndn(phi,k,nx,sn,cn,dn) !Calculate elliptic functions for all ! values of phi tv1=lapphi - etabar*v + noise & - 2.0*n*(n+1)*k*k*sn*cn*dn tphi1=v phi=phibefore+eps*tphi1 v=v+eps*tv1 lapphi=(cshift(phi,-1)+cshift(phi,1)-2.0*phi)/length/length ! changed potential here also... call sncndn(phi,k,nx,sn,cn,dn) !Calculate elliptic functions tv2=lapphi - etabar*v + noise & - 2.0*n*(n+1)*k*k*sn*cn*dn v=pibefore+eps*0.5*(tv1+tv2) phi=phibefore+eps*0.5*(tphi1+v) end SUBROUTINE rk2 !-------SUBROUTINE sncndn------------------------------------------------------ ! sn,cn, and dn are jacobian elliptic functions ! In our case, these functions are used to calculate the potential and force. ! For us, "u" will be the value of the field, phi, while "k" determines the ! shape of the potential. k=0 reduces to a sin curve, while k=1 reduces to ! the tanh curve. SUBROUTINE sncndn(u,k,nx,sn,cn,dn) ! Jacobi elliptic functions by AGM method implicit none ! Dummy variables integer, intent(in):: nx real, dimension(nx), intent(in):: u real, dimension(nx), intent(out):: sn,cn,dn !Contains the results. real, intent(in):: k !hpf$ align (:) with u(:) :: sn,cn,dn ! Local variables integer :: i !Loop variable real :: eps=10E-8 !Determines the precision real, dimension(0:9):: a,b,c ! we need a 9 element array for each element of the array u, in order to ! calculate an intermediate value. We create a 2 dimensional array. real, dimension(0:9,nx):: eeks !hpf& align eeks(*,:) with u(:) integer :: n,max=10 a(0)=1.0 b(0)=sqrt(1.0-k*k) c(0)=sqrt(1.0*k*k) n=0 do if(.NOT.(abs(c(n)) > eps)) exit a(n+1) = 0.5*(a(n)+b(n)) b(n+1) = sqrt(a(n)*b(n)) c(n+1) = 0.5*(a(n)-b(n)) if(n >= max-2) then write(*,*) "!!In JacobiSn n=", n, " too large" write(*,*) "Increase max or decrease eps" c(max-1) = 0.0 endif n=n+1 end do n=n-1 forall(i=1:nx)eeks(n,i)=(2**n)*a(n)*u(i) do while (n>0) forall(i=1:nx)eeks(n-1,i) = 0.5*(eeks(n,i)+asin(c(n)* & sin(eeks(n,i))/a(n))) n=n-1 end do sn(:) = sin(eeks(0,:)) cn(:) = cos(eeks(0,:)) dn(:) = cn/cos(eeks(1,:)-eeks(0,:)) end SUBROUTINE sncndn !-------SUBROUTINE PERPDF------------------------------------------------------ ! The purpose of this subroutine is to "periodify" the numerical pdf. SUBROUTINE perpdf(npdf,pdfsym,pdf_sum,spacing,pdfint,naves) implicit none ! Dummy variables: ! "naves" is number of averages real,intent(in) :: spacing integer,intent(in) :: naves,npdf,pdfint real,dimension(npdf),intent(inout) :: pdf_sum real,dimension(npdf), intent(in) :: pdfsym !hpf$ distribute pdfsym(block) !hpf$ align (:) with pdfsym(:) :: pdf_sum ! Local variables: integer :: k,lspacing lspacing=0 print *, "spacing is", spacing print *, "pdfint is", pdfint ! The scaling between phi and the npdf sized lattice is given by pdfint lspacing=int(spacing*pdfint) !cast it to an integer ! check print *, "lspacing is ", lspacing !question: do I have to perform an explicit summing over the indices of the !array? print *, "performing the periodification process promptly" do k=1,naves pdf_sum=pdf_sum+ $ eoshift(pdfsym,((-((naves-1)/2)-1)+k)*lspacing) enddo ! Above, we started the averaging with the left-most shifted pdf, then ! worked our way right until adding all of the pdfs ! Perform the average over the number of images summed up pdf_sum=pdf_sum/naves end SUBROUTINE perpdf