diff --git a/Code/Makefile.gfortran.openmp b/Code/Makefile.gfortran.openmp
new file mode 100644
index 0000000000000000000000000000000000000000..e69bdd1bc7258e584551b8476c03b65768a12e11
--- /dev/null
+++ b/Code/Makefile.gfortran.openmp
@@ -0,0 +1,78 @@
+COMPILER = gfortran
+COMPILERFLAGS = -O3 -fopenmp -msse4.2 -mfpmath=sse -ffast-math -finline-functions -funroll-loops
+#COMPILERFLAGS = -g -fbacktrace
+LINKERFLAGS = $(COMPILERFLAGS)
+
+SOURCE_TRAILER = .f90
+EXEC = sky3d.omp
+OBJS = params.o grids.o levels.o fourier.o forces.o sequential.o \
+ inout.o coulomb.o trivial.o densities.o fragments.o twobody.o \
+ energies.o static.o meanfield.o dynamic.o pairs.o moment.o \
+ main3d.o user.o external.o
+
+SHELL = /bin/sh
+
+#---------------------------------------------------------
+all:$(EXEC)
+$(EXEC):$(OBJS)
+ $(COMPILER) $(LINKERFLAGS) -o $@ $(OBJS) -lfftw3 -llapack
+#
+
+user.o : user.f90 params.o grids.o levels.o
+
+external.o : external.f90 params.o sequential.o grids.o levels.o meanfield.o
+
+coulomb.o : coulomb.f90 params.o grids.o densities.o
+
+inout.o : inout.f90 params.o sequential.o grids.o moment.o densities.o \
+ levels.o coulomb.o forces.o meanfield.o forces.o
+
+trivial.o : trivial.f90 params.o grids.o
+
+forces.o : forces.f90 forces.data params.o
+
+fourier.o : fourier.f90 params.o grids.o
+
+params.o : params.f90
+
+moment.o : moment.f90 params.o grids.o densities.o
+
+sequential.o: sequential.f90 levels.o params.o grids.o densities.o
+
+levels.o : levels.f90 params.o grids.o forces.o fourier.o trivial.o
+
+grids.o : grids.f90 params.o forces.o
+
+densities.o : densities.f90 params.o grids.o levels.o trivial.o
+
+fragments.o : fragments.f90 params.o grids.o forces.o levels.o twobody.o \
+ sequential.o
+
+twobody.o : twobody.f90 params.o grids.o densities.o moment.o forces.o
+
+energies.o : energies.f90 params.o forces.o densities.o levels.o grids.o \
+ trivial.o pairs.o
+
+meanfield.o : meanfield.f90 params.o densities.o forces.o grids.o coulomb.o \
+ trivial.o levels.o
+
+pairs.o : pairs.f90 forces.o grids.o levels.o densities.o
+
+static.o : static.f90 params.o densities.o grids.o \
+ levels.o moment.o energies.o inout.o pairs.o meanfield.o
+
+dynamic.o : dynamic.f90 params.o densities.o grids.o trivial.o \
+ levels.o moment.o energies.o inout.o meanfield.o sequential.o \
+ twobody.o external.o
+
+main3d.o : main3d.f90 params.o fourier.o forces.o densities.o meanfield.o levels.o grids.o fragments.o sequential.o dynamic.o static.o coulomb.o user.o
+
+.SUFFIXES: .f90 .f .o
+
+%.o : %.mod
+
+.f90.o:
+ $(COMPILER) $(COMPILERFLAGS) -c $<
+
+.f.o:
+ $(COMPILER) $(COMPILERFLAGS) -c $<
diff --git a/Code/coulomb.f90 b/Code/coulomb.f90
index 1f97a7941289f400a32be2ad47070e826ad2aca8..b56fbd06eb6baf8a6eb08e4367edb6ebaaecfdb5 100644
--- a/Code/coulomb.f90
+++ b/Code/coulomb.f90
@@ -11,7 +11,7 @@ MODULE Coulomb
PUBLIC :: poisson,coulinit,wcoul
PRIVATE :: initiq
CONTAINS
-!***************************************************
+ !***************************************************
SUBROUTINE poisson
COMPLEX(db) :: rho2(nx2,ny2,nz2)
! put proton density into array of same or double size, zeroing rest
@@ -30,8 +30,8 @@ CONTAINS
CALL dfftw_execute_dft(coulplan2,rho2,rho2)
wcoul=REAL(rho2(1:nx,1:ny,1:nz))/(nx2*ny2*nz2)
END SUBROUTINE poisson
-!***************************************************
- SUBROUTINE coulinit
+ !***************************************************
+ SUBROUTINE coulinit
INCLUDE 'fftw3.f'
REAL(db),ALLOCATABLE :: iqx(:),iqy(:),iqz(:)
INTEGER :: i,j,k
@@ -66,7 +66,7 @@ CONTAINS
END IF
DEALLOCATE(iqx,iqy,iqz)
END SUBROUTINE coulinit
-!***************************************************
+ !***************************************************
SUBROUTINE initiq(n,d,iq)
INTEGER,INTENT(IN) :: n
REAL(db),INTENT(IN) :: d
diff --git a/Code/densities.f90 b/Code/densities.f90
index 210b7a12c0dcb2a72cce2e0ac453e9885df1311c..4446bc4da5dfe2301c8a143d98b7309ff6785bcc 100644
--- a/Code/densities.f90
+++ b/Code/densities.f90
@@ -30,21 +30,21 @@ CONTAINS
(psin(:,:,:,1)*CONJG(psin(:,:,:,1))+ &
psin(:,:,:,2)*CONJG(psin(:,:,:,2)))
FORALL(iz=1:nz,iy=1:ny,ix=1:nx)
- lsdens(ix,iy,iz,1,iq)=lsdens(ix,iy,iz,1,iq)+2.D0*weight &
- *REAL(CONJG(psin(ix,iy,iz,1))*psin(ix,iy,iz,2))
- lsdens(ix,iy,iz,2,iq)=lsdens(ix,iy,iz,2,iq)+2.D0*weight &
- *AIMAG(CONJG(psin(ix,iy,iz,1))*psin(ix,iy,iz,2))
- lsdens(ix,iy,iz,3,iq)=lsdens(ix,iy,iz,3,iq)+weight &
- *(REAL(CONJG(psin(ix,iy,iz,1))*psin(ix,iy,iz,1)) &
- -REAL(CONJG(psin(ix,iy,iz,2))*psin(ix,iy,iz,2)))
+ lsdens(ix,iy,iz,1,iq)=lsdens(ix,iy,iz,1,iq)+2.D0*weight &
+ *REAL(CONJG(psin(ix,iy,iz,1))*psin(ix,iy,iz,2))
+ lsdens(ix,iy,iz,2,iq)=lsdens(ix,iy,iz,2,iq)+2.D0*weight &
+ *AIMAG(CONJG(psin(ix,iy,iz,1))*psin(ix,iy,iz,2))
+ lsdens(ix,iy,iz,3,iq)=lsdens(ix,iy,iz,3,iq)+weight &
+ *(REAL(CONJG(psin(ix,iy,iz,1))*psin(ix,iy,iz,1)) &
+ -REAL(CONJG(psin(ix,iy,iz,2))*psin(ix,iy,iz,2)))
END FORALL
!***********************************************************************
! x-derivatives
!***********************************************************************
IF(TFFT) THEN
- CALL cdervx(psin,ps1)
+ CALL cdervx(psin,ps1)
ELSE
- CALL cmulx(der1x,psin,ps1,0)
+ CALL cmulx(der1x,psin,ps1,0)
ENDIF
ltau(:,:,:,iq)=ltau(:,:,:,iq)+weight* &
(ps1(:,:,:,1)*CONJG(ps1(:,:,:,1))+ps1(:,:,:,2)*CONJG(ps1(:,:,:,2)))
@@ -52,20 +52,20 @@ CONTAINS
AIMAG(ps1(:,:,:,1)*CONJG(psin(:,:,:,1))+ &
ps1(:,:,:,2)*CONJG(psin(:,:,:,2)))
FORALL(ix=1:nx,iy=1:ny,iz=1:nz)
- lsodens(ix,iy,iz,2,iq)=lsodens(ix,iy,iz,2,iq)-weight &
- *( AIMAG(ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,1))) &
- - AIMAG(ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,2))) )
- lsodens(ix,iy,iz,3,iq)=lsodens(ix,iy,iz,3,iq)-weight &
- *( REAL(psin(ix,iy,iz,1)*CONJG(ps1(ix,iy,iz,2))) &
- - REAL(psin(ix,iy,iz,2)*CONJG(ps1(ix,iy,iz,1))) )
+ lsodens(ix,iy,iz,2,iq)=lsodens(ix,iy,iz,2,iq)-weight &
+ *( AIMAG(ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,1))) &
+ - AIMAG(ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,2))) )
+ lsodens(ix,iy,iz,3,iq)=lsodens(ix,iy,iz,3,iq)-weight &
+ *( REAL(psin(ix,iy,iz,1)*CONJG(ps1(ix,iy,iz,2))) &
+ - REAL(psin(ix,iy,iz,2)*CONJG(ps1(ix,iy,iz,1))) )
END FORALL
!***********************************************************************
! y-derivatives
!***********************************************************************
IF(TFFT) THEN
- CALL cdervy(psin,ps1)
+ CALL cdervy(psin,ps1)
ELSE
- CALL cmuly(der1y,psin,ps1,0)
+ CALL cmuly(der1y,psin,ps1,0)
ENDIF
ltau(:,:,:,iq)=ltau(:,:,:,iq)+weight* &
(ps1(:,:,:,1)*CONJG(ps1(:,:,:,1))+ps1(:,:,:,2)*CONJG(ps1(:,:,:,2)))
@@ -73,20 +73,20 @@ CONTAINS
AIMAG(ps1(:,:,:,1)*CONJG(psin(:,:,:,1))+ &
ps1(:,:,:,2)*CONJG(psin(:,:,:,2)))
FORALL(ix=1:Nx,iy=1:ny,iz=1:nz)
- lsodens(ix,iy,iz,1,iq)=lsodens(ix,iy,iz,1,iq)+weight &
- *AIMAG( ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,1)) &
- -ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,2)) )
- lsodens(ix,iy,iz,3,iq)=lsodens(ix,iy,iz,3,iq)-weight &
- *AIMAG( ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,1)) &
- +ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,2)) )
+ lsodens(ix,iy,iz,1,iq)=lsodens(ix,iy,iz,1,iq)+weight &
+ *AIMAG( ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,1)) &
+ -ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,2)) )
+ lsodens(ix,iy,iz,3,iq)=lsodens(ix,iy,iz,3,iq)-weight &
+ *AIMAG( ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,1)) &
+ +ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,2)) )
END FORALL
!***********************************************************************
! z-derivatives
!***********************************************************************
IF(TFFT) THEN
- CALL cdervz(psin,ps1)
+ CALL cdervz(psin,ps1)
ELSE
- CALL cmulz(der1z,psin,ps1,0)
+ CALL cmulz(der1z,psin,ps1,0)
ENDIF
ltau(:,:,:,iq)=ltau(:,:,:,iq)+weight* &
(ps1(:,:,:,1)*CONJG(ps1(:,:,:,1))+ps1(:,:,:,2)*CONJG(ps1(:,:,:,2)))
@@ -94,12 +94,12 @@ CONTAINS
AIMAG(ps1(:,:,:,1)*CONJG(psin(:,:,:,1))+ &
ps1(:,:,:,2)*CONJG(psin(:,:,:,2)))
FORALL(ix=1:Nx,iy=1:ny,iz=1:nz)
- lsodens(ix,iy,iz,1,iq)=lsodens(ix,iy,iz,1,iq)+weight &
- *REAL( ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,1)) &
- -ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,2)) )
- lsodens(ix,iy,iz,2,iq)=lsodens(ix,iy,iz,2,iq)+weight &
- *AIMAG( ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,1)) &
- +ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,2)) )
+ lsodens(ix,iy,iz,1,iq)=lsodens(ix,iy,iz,1,iq)+weight &
+ *REAL( ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,1)) &
+ -ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,2)) )
+ lsodens(ix,iy,iz,2,iq)=lsodens(ix,iy,iz,2,iq)+weight &
+ *AIMAG( ps1(ix,iy,iz,2)*CONJG(psin(ix,iy,iz,1)) &
+ +ps1(ix,iy,iz,1)*CONJG(psin(ix,iy,iz,2)) )
END FORALL
END SUBROUTINE add_density
END MODULE Densities
diff --git a/Code/dynamic.f90 b/Code/dynamic.f90
index 66984ac85d6b7495da2981453a58c994d9759a6d..f7b6208dfb8d03b6b5e51a97bea5b72b8fe59a28 100644
--- a/Code/dynamic.f90
+++ b/Code/dynamic.f90
@@ -28,12 +28,12 @@ CONTAINS
NAMELIST /dynamic/ nt,dt,mxpact,mrescm,rsep,texternal
READ(5,dynamic)
IF(wflag) THEN
- WRITE(*,*) '***** Parameters for the dynamic calculation *****'
- WRITE(*,"(A,I8,A,F10.5,A)") " Number of time steps:",nt, &
- ' Time step size: ',dt,' fm/c'
- WRITE(*,'(A,F7.2,A)') ' The calculation stops at ',rsep, &
- ' fm fragment separation'
- WRITE(*,'(A,I3)') ' Power limit in operator expansion:',mxpact
+ WRITE(*,*) '***** Parameters for the dynamic calculation *****'
+ WRITE(*,"(A,I8,A,F10.5,A)") " Number of time steps:",nt, &
+ ' Time step size: ',dt,' fm/c'
+ WRITE(*,'(A,F7.2,A)') ' The calculation stops at ',rsep, &
+ ' fm fragment separation'
+ WRITE(*,'(A,I3)') ' Power limit in operator expansion:',mxpact
ENDIF
IF(texternal) CALL getin_external
END SUBROUTINE getin_dynamic
@@ -43,37 +43,37 @@ CONTAINS
COMPLEX(db) :: ps4(nx,ny,nz,2)
! Step 1: Preparation phase
IF(.NOT.trestart) THEN
- iter=0
- time=0.0D0
- ! save wave functions
- CALL write_wavefunctions
- IF(wflag) WRITE(*,*) 'Wrote wave function file after initialization'
+ iter=0
+ time=0.0D0
+ ! save wave functions
+ CALL write_wavefunctions
+ IF(wflag) WRITE(*,*) 'Wrote wave function file after initialization'
END IF
! external boost
IF(texternal) THEN
- CALL extboost(text_timedep)
+ CALL extboost(text_timedep)
END IF
! Create protocol files
IF(wflag) THEN
- ! Initialize *.res files
- CALL start_protocol(energiesfile, &
- '# Time N(n) N(p) E(sum) E(integ) Ekin &
- & Ecoll(n) Ecoll(p)')
- CALL start_protocol(monopolesfile, &
- '# Time rms_n rms_p rms_tot rms_n-rms_p')
- CALL start_protocol(dipolesfile, &
- '# Iter c.m. x-y-z Isovector&
- &dipoles x-y-z')
- CALL start_protocol(quadrupolesfile, &
- '# Time Q(n) Q(p) Q(n+p) x²(n) &
- &y²(n) z²(n) x²(p) y²(p) z²(p)')
- CALL start_protocol(spinfile, &
- '# Iter Lx Ly Lz Sx Sy &
- &Sz Jx Jy Jz')
- CALL start_protocol(momentafile, &
- '# Time Px Py Pz')
- IF(texternal) CALL start_protocol(extfieldfile, &
- '# time average_extfield')
+ ! Initialize *.res files
+ CALL start_protocol(energiesfile, &
+ '# Time N(n) N(p) E(sum) E(integ) Ekin &
+ & Ecoll(n) Ecoll(p)')
+ CALL start_protocol(monopolesfile, &
+ '# Time rms_n rms_p rms_tot rms_n-rms_p')
+ CALL start_protocol(dipolesfile, &
+ '# Iter c.m. x-y-z Isovector&
+ &dipoles x-y-z')
+ CALL start_protocol(quadrupolesfile, &
+ '# Time Q(n) Q(p) Q(n+p) x²(n) &
+ &y²(n) z²(n) x²(p) y²(p) z²(p)')
+ CALL start_protocol(spinfile, &
+ '# Iter Lx Ly Lz Sx Sy &
+ &Sz Jx Jy Jz')
+ CALL start_protocol(momentafile, &
+ '# Time Px Py Pz')
+ IF(texternal) CALL start_protocol(extfieldfile, &
+ '# time average_extfield')
END IF
! calculate densities and currents
@@ -85,8 +85,8 @@ CONTAINS
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) &
!$OMP REDUCTION(+:rho,tau,current,sdens,sodens)
DO nst=1,nstloc
- CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), &
- psi(:,:,:,:,nst),rho,tau,current,sdens,sodens)
+ CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), &
+ psi(:,:,:,:,nst),rho,tau,current,sdens,sodens)
ENDDO
!$OMP END PARALLEL DO
IF(tmpi) CALL collect_densities
@@ -98,82 +98,82 @@ CONTAINS
istart=iter+1
! Step 2: start loop and do half-time step
Timestepping: DO iter=istart,nt
- IF(wflag) WRITE(*,'(/A,I6,A,F8.2,A)') ' Starting time step #',iter, &
- ' at time=',time,' fm/c'
- ! correction for parallel version
- IF(tmpi) THEN
- rho=rho/mpi_nprocs
- tau=tau/mpi_nprocs
- current=current/mpi_nprocs
- sodens=sodens/mpi_nprocs
- sdens=sdens/mpi_nprocs
- ENDIF
- ! propagate to end of time step and add to densities
- !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,ps4) SCHEDULE(STATIC) &
- !$OMP REDUCTION(+:rho,tau,current,sdens,sodens)
- DO nst=1,nstloc
- ps4=psi(:,:,:,:,nst)
- CALL tstep(isospin(globalindex(nst)),mxpact/2,ps4)
- CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), &
- ps4,rho,tau,current,sdens,sodens)
- ENDDO
- !$OMP END PARALLEL DO
- IF(tmpi) CALL collect_densities
- ! average over time step
- rho=0.5D0*rho
- tau=0.5D0*tau
- current=0.5D0*current
- sodens=0.5D0*sodens
- sdens=0.5D0*sdens
- ! compute mean field and add external field
- CALL skyrme
- IF(text_timedep) CALL extfld(time+dt/2.0D0)
- ! Step 3: full time step
- ! reset densities
- rho=0.0D0
- tau=0.0D0
- current=0.0D0
- sdens=0.0D0
- sodens=0.0D0
- ! propagate to end of step, accumulate densities
- !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,ps4) SCHEDULE(STATIC) &
- !$OMP REDUCTION(+:rho,tau,current,sdens,sodens)
- DO nst=1,nstloc
- ps4=psi(:,:,:,:,nst)
- CALL tstep(isospin(globalindex(nst)),mxpact,ps4)
- CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), &
- ps4,rho,tau,current,sdens,sodens)
- psi(:,:,:,:,nst)=ps4
- ENDDO
- !$OMP END PARALLEL DO
- ! sum up over nodes
- IF(tmpi) CALL collect_densities
- ! Step 4: eliminate center-of-mass motion if desired
- IF(mrescm/=0) THEN
- IF(MOD(iter,mrescm)==0) THEN
- CALL resetcm
- !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) &
- !$OMP REDUCTION(+:rho,tau,current,sdens,sodens)
- DO nst=1,nstloc
- CALL add_density(isospin(globalindex(nst)), &
- wocc(globalindex(nst)), &
- psi(:,:,:,:,nst),rho,tau,current,sdens,sodens)
- ENDDO
- !$OMP END PARALLEL DO
- IF(tmpi) CALL collect_densities
- ENDIF
+ IF(wflag) WRITE(*,'(/A,I6,A,F8.2,A)') ' Starting time step #',iter, &
+ ' at time=',time,' fm/c'
+ ! correction for parallel version
+ IF(tmpi) THEN
+ rho=rho/mpi_nprocs
+ tau=tau/mpi_nprocs
+ current=current/mpi_nprocs
+ sodens=sodens/mpi_nprocs
+ sdens=sdens/mpi_nprocs
+ ENDIF
+ ! propagate to end of time step and add to densities
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,ps4) SCHEDULE(STATIC) &
+ !$OMP REDUCTION(+:rho,tau,current,sdens,sodens)
+ DO nst=1,nstloc
+ ps4=psi(:,:,:,:,nst)
+ CALL tstep(isospin(globalindex(nst)),mxpact/2,ps4)
+ CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), &
+ ps4,rho,tau,current,sdens,sodens)
+ ENDDO
+ !$OMP END PARALLEL DO
+ IF(tmpi) CALL collect_densities
+ ! average over time step
+ rho=0.5D0*rho
+ tau=0.5D0*tau
+ current=0.5D0*current
+ sodens=0.5D0*sodens
+ sdens=0.5D0*sdens
+ ! compute mean field and add external field
+ CALL skyrme
+ IF(text_timedep) CALL extfld(time+dt/2.0D0)
+ ! Step 3: full time step
+ ! reset densities
+ rho=0.0D0
+ tau=0.0D0
+ current=0.0D0
+ sdens=0.0D0
+ sodens=0.0D0
+ ! propagate to end of step, accumulate densities
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,ps4) SCHEDULE(STATIC) &
+ !$OMP REDUCTION(+:rho,tau,current,sdens,sodens)
+ DO nst=1,nstloc
+ ps4=psi(:,:,:,:,nst)
+ CALL tstep(isospin(globalindex(nst)),mxpact,ps4)
+ CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), &
+ ps4,rho,tau,current,sdens,sodens)
+ psi(:,:,:,:,nst)=ps4
+ ENDDO
+ !$OMP END PARALLEL DO
+ ! sum up over nodes
+ IF(tmpi) CALL collect_densities
+ ! Step 4: eliminate center-of-mass motion if desired
+ IF(mrescm/=0) THEN
+ IF(MOD(iter,mrescm)==0) THEN
+ CALL resetcm
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) &
+ !$OMP REDUCTION(+:rho,tau,current,sdens,sodens)
+ DO nst=1,nstloc
+ CALL add_density(isospin(globalindex(nst)), &
+ wocc(globalindex(nst)), &
+ psi(:,:,:,:,nst),rho,tau,current,sdens,sodens)
+ ENDDO
+ !$OMP END PARALLEL DO
+ IF(tmpi) CALL collect_densities
+ ENDIF
+ ENDIF
+ ! Step 5: generating some output
+ time=time+dt
+ CALL tinfo
+ ! Step 6: finishing up
+ ! compute densities, currents, potentials etc. *
+ CALL skyrme
+ IF(text_timedep) CALL extfld(time+dt)
+ IF(MOD(iter,mrest)==0) THEN
+ CALL write_wavefunctions
+ IF(wflag) WRITE(*,*) ' Wrote restart file at end of iter=',iter
ENDIF
- ! Step 5: generating some output
- time=time+dt
- CALL tinfo
- ! Step 6: finishing up
- ! compute densities, currents, potentials etc. *
- CALL skyrme
- IF(text_timedep) CALL extfld(time+dt)
- IF(MOD(iter,mrest)==0) THEN
- CALL write_wavefunctions
- IF(wflag) WRITE(*,*) ' Wrote restart file at end of iter=',iter
- ENDIF
END DO Timestepping
END SUBROUTINE dynamichf
!***********************************************************************
@@ -188,20 +188,20 @@ CONTAINS
! compute exp(-I*dt*h) by power series expansion *
!***********************************************************************
DO m=1,mxp
- fmd=-dt/(hbc*m)
- CALL hpsi(iq,esf,ps1,ps2)
- DO is=1,2
- DO iz=1,nz
- DO iy=1,ny
- DO ix=1,nx
- ps1(ix,iy,iz,is)=&
- CMPLX(-fmd*AIMAG(ps2(ix,iy,iz,is)), &
- fmd*REAL(ps2(ix,iy,iz,is)),db)
- psout(ix,iy,iz,is)=psout(ix,iy,iz,is)+ps1(ix,iy,iz,is)
- ENDDO
+ fmd=-dt/(hbc*m)
+ CALL hpsi(iq,esf,ps1,ps2)
+ DO is=1,2
+ DO iz=1,nz
+ DO iy=1,ny
+ DO ix=1,nx
+ ps1(ix,iy,iz,is)=&
+ CMPLX(-fmd*AIMAG(ps2(ix,iy,iz,is)), &
+ fmd*REAL(ps2(ix,iy,iz,is)),db)
+ psout(ix,iy,iz,is)=psout(ix,iy,iz,is)+ps1(ix,iy,iz,is)
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
+ ENDDO
END DO
END SUBROUTINE tstep
!***********************************************************************
@@ -252,9 +252,9 @@ CONTAINS
IF(tmpi) THEN
CALL collect_sp_properties
ENDIF
- END IF
-! Step 5: total energies & angular momenta
- IF(printnow.AND.wflag) THEN
+ END IF
+ ! Step 5: total energies & angular momenta
+ IF(printnow.AND.wflag) THEN
CALL integ_energy
CALL sum_energy
DO iq=1,2
@@ -303,10 +303,10 @@ CONTAINS
ENDIF
! Step 8: check whether final distance is reached in twobody case
IF(istwobody.AND.roft>rsep.AND.rdot>=0.D0) THEN
- CALL twobody_print
- CALL write_wavefunctions
- IF(wflag) WRITE(*,*) ' Final separation distance reached'
- STOP ' Final distance reached'
+ CALL twobody_print
+ CALL write_wavefunctions
+ IF(wflag) WRITE(*,*) ' Final separation distance reached'
+ STOP ' Final distance reached'
ENDIF
initialcall=.FALSE.
END SUBROUTINE tinfo
diff --git a/Code/energies.f90 b/Code/energies.f90
index 5d2b817803d7cc48be8a6543f1ca7f2d8d1bd02b..a045175f6342d0103b6fdae1212e8c50b93e14b3 100644
--- a/Code/energies.f90
+++ b/Code/energies.f90
@@ -24,7 +24,7 @@ MODULE Energies
REAL(db) :: e3corr ! rearrangement energy
REAL(db) :: orbital(3),spin(3),total_angmom(3)
CONTAINS
-!***************************************************************************
+ !***************************************************************************
SUBROUTINE integ_energy
USE Trivial, ONLY: rmulx,rmuly,rmulz
USE Grids, ONLY: wxyz,der1x,der2x,der1y,der2y,der1z,der2z
@@ -35,29 +35,29 @@ CONTAINS
REAL(db) :: workb(nx,ny,nz,3,2)
! Step 1: compute laplacian of densities, then ehf0, ehf2, and ehf3
DO iq=1,2
- CALL rmulx(der2x,rho(:,:,:,iq),worka(:,:,:,iq),0)
- CALL rmuly(der2y,rho(:,:,:,iq),worka(:,:,:,iq),1)
- CALL rmulz(der2z,rho(:,:,:,iq),worka(:,:,:,iq),1)
+ CALL rmulx(der2x,rho(:,:,:,iq),worka(:,:,:,iq),0)
+ CALL rmuly(der2y,rho(:,:,:,iq),worka(:,:,:,iq),1)
+ CALL rmulz(der2z,rho(:,:,:,iq),worka(:,:,:,iq),1)
ENDDO
ehf0=0.0D0
ehf3=0.0D0
ehf2=0.0D0
DO iz=1,nz
- DO iy=1,ny
- DO ix=1,nx
- rhot=rho(ix,iy,iz,1)+rho(ix,iy,iz,2)
- rhon=rho(ix,iy,iz,1)
- rhop=rho(ix,iy,iz,2)
- ehf0=ehf0+wxyz*(b0*rhot**2-b0p*(rhop**2+rhon**2))/2.D0
- ehf3=ehf3+wxyz*rhot**f%power*(b3*rhot**2 &
- -b3p*(rhop**2+rhon**2))/3.D0
- d2rho=worka(ix,iy,iz,1)+worka(ix,iy,iz,2)
- d2rhon=worka(ix,iy,iz,1)
- d2rhop=worka(ix,iy,iz,2)
- ehf2=ehf2+wxyz*(-b2*rhot*d2rho+b2p*(rhop* &
- d2rhop+rhon*d2rhon))/2.D0
- ENDDO
- ENDDO
+ DO iy=1,ny
+ DO ix=1,nx
+ rhot=rho(ix,iy,iz,1)+rho(ix,iy,iz,2)
+ rhon=rho(ix,iy,iz,1)
+ rhop=rho(ix,iy,iz,2)
+ ehf0=ehf0+wxyz*(b0*rhot**2-b0p*(rhop**2+rhon**2))/2.D0
+ ehf3=ehf3+wxyz*rhot**f%power*(b3*rhot**2 &
+ -b3p*(rhop**2+rhon**2))/3.D0
+ d2rho=worka(ix,iy,iz,1)+worka(ix,iy,iz,2)
+ d2rhon=worka(ix,iy,iz,1)
+ d2rhop=worka(ix,iy,iz,2)
+ ehf2=ehf2+wxyz*(-b2*rhot*d2rho+b2p*(rhop* &
+ d2rhop+rhon*d2rhon))/2.D0
+ ENDDO
+ ENDDO
ENDDO
e3corr=-f%power*ehf3/2.D0
! Step 2: b1-contribution. worka=square of the current vector
@@ -71,9 +71,9 @@ CONTAINS
! Step 3: the spin-orbit contribution
! (a) Time even part: worka=div J
DO iq=1,2
- CALL rmulx(der1x,sodens(:,:,:,1,iq),worka(:,:,:,iq),0)
- CALL rmuly(der1y,sodens(:,:,:,2,iq),worka(:,:,:,iq),1)
- CALL rmulz(der1z,sodens(:,:,:,3,iq),worka(:,:,:,iq),1)
+ CALL rmulx(der1x,sodens(:,:,:,1,iq),worka(:,:,:,iq),0)
+ CALL rmuly(der1y,sodens(:,:,:,2,iq),worka(:,:,:,iq),1)
+ CALL rmulz(der1z,sodens(:,:,:,3,iq),worka(:,:,:,iq),1)
ENDDO
ehfls=wxyz*SUM(-b4*(rho(:,:,:,1)+rho(:,:,:,2)) &
*(worka(:,:,:,1)+worka(:,:,:,2)) &
@@ -98,28 +98,28 @@ CONTAINS
ehfc=0.0D0
ecorc=0.0D0
IF(tcoul) THEN
- IF(f%ex/=0) THEN
- sc=-3.0D0/4.0D0*slate
- ELSE
- sc=0.D0
- END IF
- DO iz=1,nz
- DO iy=1,ny
- DO ix=1,nx
- rhop=rho(ix,iy,iz,2)
- ehfc=ehfc+wxyz *(0.5D0*rhop*wcoul(ix,iy,iz) &
- +sc*rhop**(4.0D0/3.0D0))
- ecorc=ecorc+wxyz*sc/3.0D0*rhop**(4.0D0/3.0D0)
+ IF(f%ex/=0) THEN
+ sc=-3.0D0/4.0D0*slate
+ ELSE
+ sc=0.D0
+ END IF
+ DO iz=1,nz
+ DO iy=1,ny
+ DO ix=1,nx
+ rhop=rho(ix,iy,iz,2)
+ ehfc=ehfc+wxyz *(0.5D0*rhop*wcoul(ix,iy,iz) &
+ +sc*rhop**(4.0D0/3.0D0))
+ ecorc=ecorc+wxyz*sc/3.0D0*rhop**(4.0D0/3.0D0)
+ ENDDO
ENDDO
- ENDDO
- ENDDO
+ ENDDO
ENDIF
! Step 5: kinetic energy contribution
ehft=wxyz*SUM(f%h2m(1)*tau(:,:,:,1)+f%h2m(2)*tau(:,:,:,2))
! Step 6: form total energy
ehfint=ehft+ehf0+ehf1+ehf2+ehf3+ehfls+ehfc-epair(1)-epair(2)
END SUBROUTINE integ_energy
-!***************************************************************************
+ !***************************************************************************
SUBROUTINE sum_energy
USE Moment, ONLY: pnrtot
INTEGER :: i
@@ -129,11 +129,11 @@ CONTAINS
efluct1=SUM(wocc*sp_efluct1)/pnrtot
efluct2=SUM(wocc*sp_efluct2)/pnrtot
DO i=1,3
- orbital(i)=SUM(wocc*sp_orbital(i,:))
- spin(i)=SUM(wocc*sp_spin(i,:))
- total_angmom(i)=orbital(i)+spin(i)
+ orbital(i)=SUM(wocc*sp_orbital(i,:))
+ spin(i)=SUM(wocc*sp_spin(i,:))
+ total_angmom(i)=orbital(i)+spin(i)
END DO
END SUBROUTINE sum_energy
-!***************************************************************************
+ !***************************************************************************
END MODULE Energies
diff --git a/Code/external.f90 b/Code/external.f90
index b894c08e00ab628a24db5fa7f20445f479e8c770..3a2c8445981d996a38ea2d282280f88da011de8d 100644
--- a/Code/external.f90
+++ b/Code/external.f90
@@ -19,7 +19,7 @@ CONTAINS
CHARACTER(14),PARAMETER :: pulsetype(0:2)=(/ 'Instantaneous ', &
'Gaussian ','Cosine squared' /)
NAMELIST/extern/ amplq0,radext,widext,isoext,ipulse,omega,tau0,taut, &
- textfield_periodic
+ textfield_periodic
READ(5,extern)
IF(ipulse<0.OR.ipulse>2) STOP &
' External field: called with invalid pulse type'
@@ -57,11 +57,11 @@ CONTAINS
DO iy=1,ny
DO ix=1,nx
IF(textfield_periodic) THEN ! strictly periodic version
- facr=amplq0 *(2.D0*SIN(z(iz)*PI/zlim)**2 &
- -SIN(x(ix)*PI/xlim)**2-SIN(y(iy)*PI/ylim)**2)
+ facr=amplq0 *(2.D0*SIN(z(iz)*PI/zlim)**2 &
+ -SIN(x(ix)*PI/xlim)**2-SIN(y(iy)*PI/ylim)**2)
ELSE ! damped version
- facr=amplq0 *(2.D0*z(iz)**2-x(ix)**2-y(iy)**2) &
- /(1.0D0+EXP((SQRT(x(ix)**2+y(iy)**2+z(iz)**2)-radext)/widext))
+ facr=amplq0 *(2.D0*z(iz)**2-x(ix)**2-y(iy)**2) &
+ /(1.0D0+EXP((SQRT(x(ix)**2+y(iy)**2+z(iz)**2)-radext)/widext))
END IF
extfield(ix,iy,iz,1)=facr*facn
extfield(ix,iy,iz,2)=facr*facp
diff --git a/Code/forces.f90 b/Code/forces.f90
index 28f18bd84eaa393a81390901f30727ecea7cd51a..885d0208a69054ea4ff6ce9005c0eb3c6b4f9c08 100644
--- a/Code/forces.f90
+++ b/Code/forces.f90
@@ -50,26 +50,26 @@ CONTAINS
! seek for force in predefined ones
predefined=.FALSE.
DO i=1,nforce
- IF(TRIM(name)==TRIM(pforces(i)%name)) THEN
- predefined=.TRUE.
- f=pforces(i)
- ENDIF
+ IF(TRIM(name)==TRIM(pforces(i)%name)) THEN
+ predefined=.TRUE.
+ f=pforces(i)
+ ENDIF
END DO
IF(wflag)WRITE(*,*) '***** Force definition *****'
IF(predefined) THEN
- IF(wflag) WRITE(*,*) 'Using predefined force ',name
+ IF(wflag) WRITE(*,*) 'Using predefined force ',name
ELSEIF(h2m(1)<0.D0) THEN
- IF(wflag) WRITE(*,*) ' Force not found: ',name,' and none given in input'
- STOP
+ IF(wflag) WRITE(*,*) ' Force not found: ',name,' and none given in input'
+ STOP
ELSE
- IF(wflag) WRITE(*,*) 'Using new force ',name,' defined in input'
- IF(ABS(b4p)<1.0d-9) b4p=f%t4/2.0D0
- f%name=name
- f%ex=ex; f%zpe=zpe
- f%h2m=h2m
- f%t0=t0; f%t1=t1; f%t2=t2; f%t3=t3; f%t4=t4
- f%x0=x0; f%x1=x1;f%x2=x2;f%x3=x3;f%b4p=b4p;
- f%power=power
+ IF(wflag) WRITE(*,*) 'Using new force ',name,' defined in input'
+ IF(ABS(b4p)<1.0d-9) b4p=f%t4/2.0D0
+ f%name=name
+ f%ex=ex; f%zpe=zpe
+ f%h2m=h2m
+ f%t0=t0; f%t1=t1; f%t2=t2; f%t3=t3; f%t4=t4
+ f%x0=x0; f%x1=x1;f%x2=x2;f%x3=x3;f%b4p=b4p;
+ f%power=power
ENDIF
! turn off zpe if desired
IF(turnoff_zpe) THEN
@@ -90,36 +90,36 @@ CONTAINS
slate=(3.0D0/pi)**(1.0D0/3.0D0)*e2
! now set up pairing: first case of none
IF(TRIM(pairing)=='NONE') THEN
- p%v0prot=0.D0; p%v0neut=0.D0; p%rho0pr=0.16D0
- ipair=0
+ p%v0prot=0.D0; p%v0neut=0.D0; p%rho0pr=0.16D0
+ ipair=0
ELSE
- ! set predefined type of pairing
- IF(TRIM(pairing)=='VDI') ipair=5
- IF(TRIM(pairing)=='DDDI') ipair=6
-! IF(ipair==6) STOP 'DDDI pairing not implemented in this version'
- p%v0prot=v0prot
- p%v0neut=v0neut
- IF(ipair==6) p%rho0pr=rho0pr
- ! get predefined pairing if applicable
- IF(predefined) THEN
- IF(ipair==5) p=f%vdi
- IF(ipair==6) p=f%dddi
- ENDIF
- ! stop if this has not yielded meaningful parameters
- IF(p%v0prot<0.D0) STOP 'Pairing not defined properly'
- IF(wflag) WRITE(*,*) 'Using pairing type ',pairing,'=',ipair
+ ! set predefined type of pairing
+ IF(TRIM(pairing)=='VDI') ipair=5
+ IF(TRIM(pairing)=='DDDI') ipair=6
+ ! IF(ipair==6) STOP 'DDDI pairing not implemented in this version'
+ p%v0prot=v0prot
+ p%v0neut=v0neut
+ IF(ipair==6) p%rho0pr=rho0pr
+ ! get predefined pairing if applicable
+ IF(predefined) THEN
+ IF(ipair==5) p=f%vdi
+ IF(ipair==6) p=f%dddi
+ ENDIF
+ ! stop if this has not yielded meaningful parameters
+ IF(p%v0prot<0.D0) STOP 'Pairing not defined properly'
+ IF(wflag) WRITE(*,*) 'Using pairing type ',pairing,'=',ipair
END IF
! define average h2m and nucleon mass
h2ma=0.5D0*(f%h2m(1)+f%h2m(2))
nucleon_mass=hbc**2/(2.D0*h2ma)
IF(wflag) THEN
- WRITE(*,"(A)") " Skyrme Potential Parameters"
- WRITE(*,*) 'The force is: ',f%name
- WRITE(*,"(5(A6,F12.5))") "t0",f%t0,"t1",f%t1,"t2",f%t2,"t3",f%t3,"t4",f%t4
- WRITE(*,"(5(A6,F12.5))") "x0",f%x0,"x1",f%x1,"x2",f%x2,"x3",f%x3,"b4p",f%b4p
- WRITE(*,"(A6,F12.5)") "Power",f%power
- WRITE(*,"(A,I2)") " Pairing parameters: Option ipair:",ipair
- WRITE(*,"(3(A7,F12.5))") "v0prot",p%v0prot,"v0neut",p%v0neut,"rho0pr",p%rho0pr
+ WRITE(*,"(A)") " Skyrme Potential Parameters"
+ WRITE(*,*) 'The force is: ',f%name
+ WRITE(*,"(5(A6,F12.5))") "t0",f%t0,"t1",f%t1,"t2",f%t2,"t3",f%t3,"t4",f%t4
+ WRITE(*,"(5(A6,F12.5))") "x0",f%x0,"x1",f%x1,"x2",f%x2,"x3",f%x3,"b4p",f%b4p
+ WRITE(*,"(A6,F12.5)") "Power",f%power
+ WRITE(*,"(A,I2)") " Pairing parameters: Option ipair:",ipair
+ WRITE(*,"(3(A7,F12.5))") "v0prot",p%v0prot,"v0neut",p%v0neut,"rho0pr",p%rho0pr
ENDIF
END SUBROUTINE read_force
END MODULE Forces
diff --git a/Code/fragments.f90 b/Code/fragments.f90
index 580c939e8d11d4c7ab8ae673c3c4bd734e951c39..8f4a63cd137615539cfff742bf04381c0df0600b 100644
--- a/Code/fragments.f90
+++ b/Code/fragments.f90
@@ -28,66 +28,66 @@ CONTAINS
CHARACTER(8) :: forcename
NAMELIST /fragments/ filename,fcent,fboost,ecm,b,fix_boost
IF(trestart) THEN
- filename(1)=wffile
- fcent(:,1)=0.D0
- fboost(:,1)=0.D0
- fix_boost=.TRUE.
- IF(wflag) WRITE(*,*) 'Reading restart as one fragment from file ',wffile
+ filename(1)=wffile
+ fcent(:,1)=0.D0
+ fboost(:,1)=0.D0
+ fix_boost=.TRUE.
+ IF(wflag) WRITE(*,*) 'Reading restart as one fragment from file ',wffile
ELSE
- fix_boost=.FALSE.
- READ(5,fragments)
+ fix_boost=.FALSE.
+ READ(5,fragments)
END IF
IF(nof/=2.AND..NOT.fix_boost) THEN
- IF(wflag) WRITE(*,"(//,a,//)") "Non-fixed boost only for nof=2"
- STOP
+ IF(wflag) WRITE(*,"(//,a,//)") "Non-fixed boost only for nof=2"
+ STOP
ENDIF
IF(wflag) THEN
- IF(fix_boost) THEN
- WRITE(*,*) '***** Fixed boost values used *****'
- ELSE
- WRITE(*,*) '***** Boost values computed from twobody kinematics *****'
- END IF
+ IF(fix_boost) THEN
+ WRITE(*,*) '***** Fixed boost values used *****'
+ ELSE
+ WRITE(*,*) '***** Boost values computed from twobody kinematics *****'
+ END IF
END IF
fnewnpmin(:,1)=1
DO i=1,nof
- OPEN(UNIT=scratch,FILE=filename(i),STATUS='old',FORM=&
- 'unformatted')
- READ(scratch) iter,time,forcename,fnstmax(i),fnneut(i),fnprot(i), &
- fnumber(:,i),fnpsi(:,i),fcharge(i),fmass(i),fcmtot(:,i)
- IF(trestart) fcmtot(:,i)=0.D0
- IF(forcename/=f%name) THEN
- IF(wflag) WRITE(*,*) 'Forces do not agree: ',forcename, &
- ' : ',f%name,' Fragment:',filename(i)
- STOP
- END IF
- READ(scratch) fnx(i),fny(i),fnz(i),fdx,fdy,fdz,fwxyz
- CLOSE(UNIT=scratch)
- IF(fnx(i)>nx.OR.fny(i)>ny.OR.fnz(i)>nz) &
- CALL errormsg('Fragment dimensioning larger than present', &
- filename(i))
- IF(MAX(ABS(fdx-dx),ABS(fdy-dy),ABS(fdz-dz)) &
- >0.01D0) CALL errormsg('Grid spacing different for', &
- filename(i))
- IF(wflag) THEN
- WRITE(*,"(A,I2,2A)") " ***** Data for fragment #",i,' from file ',&
- filename(i)
- WRITE(*,"(A,1P,3E12.4)") " Location(x,y,z): ",fcent(:,i)
- WRITE(*,"(2(A,F9.4))") " Mass number: ",fmass(i), &
- ' Charge number: ',fcharge(i)
- IF(fix_boost) &
- WRITE(*,"(A,1P,3E12.4)") " Boost(x,y,z): ",fboost(:,i)
- END IF
- ! for static restart read also unoccupied wave functions
- IF(tstatic.AND.nof==1) THEN
- fnumber(1,i)=fnpsi(1,i)
- fnumber(2,i)=fnpsi(2,i)-fnpsi(1,i)
- ENDIF
- ! get isospin starting points in fragment file
- fnpmin(1,i)=1
- fnpmin(2,i)=fnpsi(1,i)+1
- IF(i<nof) &
- fnewnpmin(:,i+1)=fnewnpmin(:,i)+fnumber(:,i)
- fnewnpsi(:,i)=fnewnpmin(:,i)+fnumber(:,i)-1
+ OPEN(UNIT=scratch,FILE=filename(i),STATUS='old',FORM=&
+ 'unformatted')
+ READ(scratch) iter,time,forcename,fnstmax(i),fnneut(i),fnprot(i), &
+ fnumber(:,i),fnpsi(:,i),fcharge(i),fmass(i),fcmtot(:,i)
+ IF(trestart) fcmtot(:,i)=0.D0
+ IF(forcename/=f%name) THEN
+ IF(wflag) WRITE(*,*) 'Forces do not agree: ',forcename, &
+ ' : ',f%name,' Fragment:',filename(i)
+ STOP
+ END IF
+ READ(scratch) fnx(i),fny(i),fnz(i),fdx,fdy,fdz,fwxyz
+ CLOSE(UNIT=scratch)
+ IF(fnx(i)>nx.OR.fny(i)>ny.OR.fnz(i)>nz) &
+ CALL errormsg('Fragment dimensioning larger than present', &
+ filename(i))
+ IF(MAX(ABS(fdx-dx),ABS(fdy-dy),ABS(fdz-dz)) &
+ >0.01D0) CALL errormsg('Grid spacing different for', &
+ filename(i))
+ IF(wflag) THEN
+ WRITE(*,"(A,I2,2A)") " ***** Data for fragment #",i,' from file ',&
+ filename(i)
+ WRITE(*,"(A,1P,3E12.4)") " Location(x,y,z): ",fcent(:,i)
+ WRITE(*,"(2(A,F9.4))") " Mass number: ",fmass(i), &
+ ' Charge number: ',fcharge(i)
+ IF(fix_boost) &
+ WRITE(*,"(A,1P,3E12.4)") " Boost(x,y,z): ",fboost(:,i)
+ END IF
+ ! for static restart read also unoccupied wave functions
+ IF(tstatic.AND.nof==1) THEN
+ fnumber(1,i)=fnpsi(1,i)
+ fnumber(2,i)=fnpsi(2,i)-fnpsi(1,i)
+ ENDIF
+ ! get isospin starting points in fragment file
+ fnpmin(1,i)=1
+ fnpmin(2,i)=fnpsi(1,i)+1
+ IF(i<nof) &
+ fnewnpmin(:,i+1)=fnewnpmin(:,i)+fnumber(:,i)
+ fnewnpsi(:,i)=fnewnpmin(:,i)+fnumber(:,i)-1
END DO
! determine total wave function numbers
nneut=SUM(fnneut)
@@ -104,9 +104,9 @@ CONTAINS
mass_number=nneut+nprot
IF(.NOT.fix_boost.AND.nof==2.AND.tdynamic) CALL twobody_init
IF(nof==2) THEN
- ! Give fragment positions to 2-body analysis as initial guess
- centerx=fcent(1,1:2)
- centerz=fcent(3,1:2)
+ ! Give fragment positions to 2-body analysis as initial guess
+ centerx=fcent(1,1:2)
+ centerz=fcent(3,1:2)
END IF
END SUBROUTINE getin_fragments
!*******************************************************************
@@ -115,8 +115,8 @@ CONTAINS
IF(wflag) &
WRITE(*,*) '***** Input of fragment wave functions *****'
DO iff=1,nof
- CALL read_one_fragment(iff)
- CALL boost_fragment(iff)
+ CALL read_one_fragment(iff)
+ CALL boost_fragment(iff)
ENDDO
IF(wflag) &
WRITE(*,*) '***** All fragments loaded and boosted *****'
@@ -144,26 +144,26 @@ CONTAINS
READ(scratch) fwocc,fsp_energy,fsp_parity,fsp_norm,fsp_kinetic,fsp_efluct1
ps1=0.0D0
IF(wflag) THEN
- WRITE(*,'(A,I3)') ' ***** S.p. levels for fragment #',iff
- WRITE(*,'(A,3I5,A,2I5)') ' Iso #new #old Occup. E_sp &
- &Parity Norm E_kin E_fluct'
+ WRITE(*,'(A,I3)') ' ***** S.p. levels for fragment #',iff
+ WRITE(*,'(A,3I5,A,2I5)') ' Iso #new #old Occup. E_sp &
+ &Parity Norm E_kin E_fluct'
END IF
DO iq=1,2
- DO inew=fnewnpmin(iq,iff),fnewnpsi(iq,iff)
- iold=fnpmin(iq,iff)+inew-fnewnpmin(iq,iff)
- wocc(inew)=fwocc(iold)
- sp_energy(inew)=fsp_energy(iold)
- sp_parity(inew)=fsp_parity(iold)
- sp_norm(inew)=fsp_norm(iold)
- sp_kinetic(inew)=fsp_kinetic(iold)
- sp_efluct1(inew)=fsp_efluct1(iold)
- isospin(inew)=iq
- IF(wflag) THEN
- WRITE(*,'(3I6,6G12.5)') &
- iq,inew,iold,wocc(inew),sp_energy(inew),sp_parity(inew), &
- sp_norm(inew),sp_kinetic(inew),sp_efluct1(inew)
- ENDIF
- END DO
+ DO inew=fnewnpmin(iq,iff),fnewnpsi(iq,iff)
+ iold=fnpmin(iq,iff)+inew-fnewnpmin(iq,iff)
+ wocc(inew)=fwocc(iold)
+ sp_energy(inew)=fsp_energy(iold)
+ sp_parity(inew)=fsp_parity(iold)
+ sp_norm(inew)=fsp_norm(iold)
+ sp_kinetic(inew)=fsp_kinetic(iold)
+ sp_efluct1(inew)=fsp_efluct1(iold)
+ isospin(inew)=iq
+ IF(wflag) THEN
+ WRITE(*,'(3I6,6G12.5)') &
+ iq,inew,iold,wocc(inew),sp_energy(inew),sp_parity(inew), &
+ sp_norm(inew),sp_kinetic(inew),sp_efluct1(inew)
+ ENDIF
+ END DO
ENDDO
cmi(1)=(fcent(1,iff)-fcmtot(1,iff)-x(1)+fx(1))/(nx*dx)
cmi(2)=(fcent(2,iff)-fcmtot(2,iff)-y(1)+fy(1))/(ny*dy)
@@ -178,39 +178,39 @@ CONTAINS
multifile=ANY(fnode/=0)
! start reading wave functions
DO iq=1,2
- ! skip unoccupied neutron states in one-file case
- IF(iq==2.AND..NOT.multifile) THEN
- DO nst=fnumber(1,iff)+1,fnpsi(1,iff)
- READ(scratch)
- END DO
- END IF
- ! Read occupied states
- DO nst=1,fnumber(iq,iff)
- newnst=fnewnpmin(iq,iff)+nst-1
- oldnst=fnpmin(iq,iff)+nst-1 ! includes empty states
- IF(node(newnst)==mpi_myproc) THEN
- ipn=localindex(newnst)
- ! open correct file and position for multifile case
- IF(multifile) CALL locate(fnode(oldnst),flocalindex(oldnst))
- ps1=0.D0
- READ(scratch) ps1(1:fnx(iff),1:fny(iff),1:fnz(iff),:)
- DO is=1,2
- CALL dfftw_execute_dft(pforward,ps1(:,:,:,is),ps1(:,:,:,is))
- FORALL(ix=1:nx,iy=1:ny,iz=1:nz)
- ps1(ix,iy,iz,is)=ps1(ix,iy,iz,is)*akx(ix)*aky(iy)*akz(iz) &
- /DBLE(nx*ny*nz)
- END FORALL
- CALL dfftw_execute_dft(pbackward,ps1(:,:,:,is),ps1(:,:,:,is))
- WHERE(ABS(ps1)==0.D0)
- psi(:,:,:,:,ipn)=1.0D-20
- ELSEWHERE
- psi(:,:,:,:,ipn)=ps1
- ENDWHERE
- ENDDO
- ELSE
- IF(.NOT.multifile) READ(scratch)
- ENDIF
- ENDDO
+ ! skip unoccupied neutron states in one-file case
+ IF(iq==2.AND..NOT.multifile) THEN
+ DO nst=fnumber(1,iff)+1,fnpsi(1,iff)
+ READ(scratch)
+ END DO
+ END IF
+ ! Read occupied states
+ DO nst=1,fnumber(iq,iff)
+ newnst=fnewnpmin(iq,iff)+nst-1
+ oldnst=fnpmin(iq,iff)+nst-1 ! includes empty states
+ IF(node(newnst)==mpi_myproc) THEN
+ ipn=localindex(newnst)
+ ! open correct file and position for multifile case
+ IF(multifile) CALL locate(fnode(oldnst),flocalindex(oldnst))
+ ps1=0.D0
+ READ(scratch) ps1(1:fnx(iff),1:fny(iff),1:fnz(iff),:)
+ DO is=1,2
+ CALL dfftw_execute_dft(pforward,ps1(:,:,:,is),ps1(:,:,:,is))
+ FORALL(ix=1:nx,iy=1:ny,iz=1:nz)
+ ps1(ix,iy,iz,is)=ps1(ix,iy,iz,is)*akx(ix)*aky(iy)*akz(iz) &
+ /DBLE(nx*ny*nz)
+ END FORALL
+ CALL dfftw_execute_dft(pbackward,ps1(:,:,:,is),ps1(:,:,:,is))
+ WHERE(ABS(ps1)==0.D0)
+ psi(:,:,:,:,ipn)=1.0D-20
+ ELSEWHERE
+ psi(:,:,:,:,ipn)=ps1
+ ENDWHERE
+ ENDDO
+ ELSE
+ IF(.NOT.multifile) READ(scratch)
+ ENDIF
+ ENDDO
ENDDO
CLOSE(unit=scratch)
END SUBROUTINE read_one_fragment
@@ -222,12 +222,12 @@ CONTAINS
INTEGER :: i
CHARACTER(120) :: rsfp
IF(fileno/=presentfile) THEN
- CLOSE(scratch)
- WRITE(rsfp,'(I3.3,''.'',A)') fileno,wffile
- OPEN(scratch,FORM='unformatted',FILE=rsfp,STATUS='OLD')
- DO i=1,pos-1
- READ(scratch)
- END DO
+ CLOSE(scratch)
+ WRITE(rsfp,'(I3.3,''.'',A)') fileno,wffile
+ OPEN(scratch,FORM='unformatted',FILE=rsfp,STATUS='OLD')
+ DO i=1,pos-1
+ READ(scratch)
+ END DO
END IF
END SUBROUTINE locate
!*******************************************************************
@@ -272,21 +272,21 @@ CONTAINS
fboost(1,2)=v2*(dix*cost-diz*sint)
fboost(3,2)=v2*(dix*sint+diz*cost)
IF(wflag) THEN
- WRITE(*,*) '***** Two-body initialization *****'
- WRITE(*,"(2(A,F12.4),A)") " c. m. Energy:",ecm, &
- ' MeV. Impact parameter b: ',b,' fm'
- WRITE(*,"(A,F12.4)") " xli/hbar :",xli
- WRITE(*,*) 'Computed boost velocities in units of c'
- WRITE(*,'(A,I2,3G15.6)') (' Fragment #',i,fboost(:,i),i=1,2)
+ WRITE(*,*) '***** Two-body initialization *****'
+ WRITE(*,"(2(A,F12.4),A)") " c. m. Energy:",ecm, &
+ ' MeV. Impact parameter b: ',b,' fm'
+ WRITE(*,"(A,F12.4)") " xli/hbar :",xli
+ WRITE(*,*) 'Computed boost velocities in units of c'
+ WRITE(*,'(A,I2,3G15.6)') (' Fragment #',i,fboost(:,i),i=1,2)
END IF
!convert into kinetic energy for that direction
FORALL(i=1:2)
- fboost(:,i)=SIGN(0.5D0*fmass(i)*fboost(:,i)**2*nucleon_mass, &
- fboost(:,i))
+ fboost(:,i)=SIGN(0.5D0*fmass(i)*fboost(:,i)**2*nucleon_mass, &
+ fboost(:,i))
END FORALL
IF(wflag) THEN
- WRITE(*,*) 'Total translational kinetic energies in units of MeV'
- WRITE(*,'(A,I2,3G15.6)') (' Fragment #',i,fboost(:,i),i=1,2)
+ WRITE(*,*) 'Total translational kinetic energies in units of MeV'
+ WRITE(*,'(A,I2,3G15.6)') (' Fragment #',i,fboost(:,i),i=1,2)
END IF
END SUBROUTINE twobody_init
!*******************************************************************
diff --git a/Code/grids.f90 b/Code/grids.f90
index 87d28660d8e6fa81b181a4a70bcb20943811191c..e834b62942be71a6e462d26f6c2b614be34e9691 100644
--- a/Code/grids.f90
+++ b/Code/grids.f90
@@ -21,21 +21,21 @@ CONTAINS
dz=0.D0
READ(5,Grid)
IF(MOD(nx,2)/=0.OR.MOD(ny,2)/=0.OR.MOD(nz,2)/=0) THEN
- IF(wflag) WRITE(*,'(A,3I4)') 'Dimensions must be even: ',nx,ny,nz
- STOP
+ IF(wflag) WRITE(*,'(A,3I4)') 'Dimensions must be even: ',nx,ny,nz
+ STOP
END IF
IF(wflag) THEN
- WRITE(*,*) '***** Grid Definition *****'
- IF(periodic) THEN
- WRITE(*,*) 'Grid is periodic'
- ELSE
- WRITE(*,*) 'Grid is not periodic'
- END IF
+ WRITE(*,*) '***** Grid Definition *****'
+ IF(periodic) THEN
+ WRITE(*,*) 'Grid is periodic'
+ ELSE
+ WRITE(*,*) 'Grid is not periodic'
+ END IF
END IF
IF(dx*dy*dz<=0.D0) THEN
- IF(dx<=0.D0) STOP 'Grid spacing given as zero'
- dy=dx
- dz=dx
+ IF(dx<=0.D0) STOP 'Grid spacing given as zero'
+ dy=dx
+ dz=dx
END IF
CALL init_coord('x',nx,dx,x,der1x,der2x,cdmpx)
CALL init_coord('y',ny,dy,y,der1y,der2y,cdmpy)
diff --git a/Code/inout.f90 b/Code/inout.f90
index b89eea84c9718629d14e13f628882ff8da4bdf0a..c1cbda285b18d52381898619840869bec18d060e 100644
--- a/Code/inout.f90
+++ b/Code/inout.f90
@@ -19,30 +19,30 @@ CONTAINS
INTEGER :: nst,iq,number(2)
CHARACTER(120) :: rsfp
IF(wflag.AND.(wffile=='none'.OR.wffile=='NONE')) THEN
- WRITE(*,*) " wffile='NONE' --> no file written "
- RETURN
+ WRITE(*,*) " wffile='NONE' --> no file written "
+ RETURN
ENDIF
! Determine number of states with non-zero occupation
DO iq=1,2
- number(iq)=COUNT(wocc(npmin(iq):npsi(iq))>0.D0)
+ number(iq)=COUNT(wocc(npmin(iq):npsi(iq))>0.D0)
END DO
IF(mpi_myproc==0) THEN
- OPEN(UNIT=scratch2,FILE=wffile,STATUS='REPLACE',FORM='UNFORMATTED')
- WRITE(scratch2) iter,time,f%name,nstmax,nneut,nprot,number,npsi, &
- charge_number,mass_number,cm
- WRITE(scratch2) nx,ny,nz,dx,dy,dz,wxyz
- WRITE(scratch2) x,y,z
- WRITE(scratch2) wocc,sp_energy,sp_parity,sp_norm,sp_kinetic, &
- sp_efluct1
- WRITE(scratch2) node,localindex
- IF(mpi_nprocs>1) CLOSE(UNIT=scratch2,STATUS='KEEP')
+ OPEN(UNIT=scratch2,FILE=wffile,STATUS='REPLACE',FORM='UNFORMATTED')
+ WRITE(scratch2) iter,time,f%name,nstmax,nneut,nprot,number,npsi, &
+ charge_number,mass_number,cm
+ WRITE(scratch2) nx,ny,nz,dx,dy,dz,wxyz
+ WRITE(scratch2) x,y,z
+ WRITE(scratch2) wocc,sp_energy,sp_parity,sp_norm,sp_kinetic, &
+ sp_efluct1
+ WRITE(scratch2) node,localindex
+ IF(mpi_nprocs>1) CLOSE(UNIT=scratch2,STATUS='KEEP')
ENDIF
IF(mpi_nprocs>1) THEN
- WRITE(rsfp,'(I3.3,''.'',A)') mpi_myproc,wffile
- OPEN(UNIT=scratch2,FILE=rsfp,STATUS='REPLACE',FORM='UNFORMATTED')
+ WRITE(rsfp,'(I3.3,''.'',A)') mpi_myproc,wffile
+ OPEN(UNIT=scratch2,FILE=rsfp,STATUS='REPLACE',FORM='UNFORMATTED')
ENDIF
DO nst=1,nstloc
- WRITE(scratch2) psi(:,:,:,:,nst)
+ WRITE(scratch2) psi(:,:,:,:,nst)
ENDDO
CLOSE(UNIT=scratch2,STATUS='KEEP')
END SUBROUTINE write_wavefunctions
@@ -58,24 +58,24 @@ CONTAINS
WRITE(scratch) iter,time,nx,ny,nz
WRITE(scratch) dx,dy,dz,wxyz,x,y,z
DO i=1,nselect
- c=writeselect(i:i)
- SELECT CASE(c)
- CASE('r','R')
- CALL write_one_density('Rho',rho)
- CASE('t','T')
- CALL write_one_density('Tau',tau)
- CASE('u','U')
- CALL write_one_density('Upot',upot)
- CASE('w','W')
- WRITE(scratch) 'Wcoul ',.FALSE.,.FALSE.
- WRITE(scratch) wcoul
- CASE('c','C')
- CALL write_vec_density('Current',current)
- CASE('s','S')
- CALL write_vec_density('Spindens',sdens)
- CASE('o','O')
- CALL write_vec_density('s-o-Dens',sodens)
- END SELECT
+ c=writeselect(i:i)
+ SELECT CASE(c)
+ CASE('r','R')
+ CALL write_one_density('Rho',rho)
+ CASE('t','T')
+ CALL write_one_density('Tau',tau)
+ CASE('u','U')
+ CALL write_one_density('Upot',upot)
+ CASE('w','W')
+ WRITE(scratch) 'Wcoul ',.FALSE.,.FALSE.
+ WRITE(scratch) wcoul
+ CASE('c','C')
+ CALL write_vec_density('Current',current)
+ CASE('s','S')
+ CALL write_vec_density('Spindens',sdens)
+ CASE('o','O')
+ CALL write_vec_density('s-o-Dens',sodens)
+ END SELECT
END DO
CLOSE(UNIT=scratch)
END SUBROUTINE write_densities
@@ -88,10 +88,10 @@ CONTAINS
stored_name=name
WRITE(scratch) stored_name,.FALSE.,write_isospin
IF(write_isospin) THEN
- WRITE(scratch) values
+ WRITE(scratch) values
ELSE
- a=values(:,:,:,1)+values(:,:,:,2)
- WRITE(scratch) a
+ a=values(:,:,:,1)+values(:,:,:,2)
+ WRITE(scratch) a
END IF
END SUBROUTINE write_one_density
!**************************************************************************
@@ -103,10 +103,10 @@ CONTAINS
stored_name=name
WRITE(scratch) stored_name,.TRUE.,write_isospin
IF(write_isospin) THEN
- WRITE(scratch) values
+ WRITE(scratch) values
ELSE
- a=values(:,:,:,:,1)+values(:,:,:,:,2)
- WRITE(scratch) a
+ a=values(:,:,:,:,1)+values(:,:,:,:,2)
+ WRITE(scratch) a
END IF
END SUBROUTINE write_vec_density
!**************************************************************************
@@ -137,28 +137,28 @@ CONTAINS
WRITE(*,'(A)') ' z '
WRITE(*,'(1X,F6.2,1X,121A)') zco(1),ibor(1:nhor)
DO i=2,nver-1
- zcu=z(nz)-(i-1)*dzp
- DO j=2,nhor-1
- xcu=x(1)+(j-1)*dxp
- jcar=1+INT(0.5D0+ &
- bplina(nx,nz,x,z,rhoplt,xcu,zcu) &
- *5.D0/density_scale)
- jcar=MIN(jcar,25)
- ifun(j)=icar(jcar)
- END DO
- IF(MOD(i-1,izsc)==0) THEN
- ntz=(i+izsc-1)/izsc
- ifun(nhor)='+'
- WRITE(*,'(1X,F6.2,1X,A,120A)') zco(ntz),'+',ifun(2:nhor)
- ELSE
- ifun(nhor)='i'
- WRITE(*,'(8X,A,120A)') 'i',ifun(2:nhor)
- END IF
+ zcu=z(nz)-(i-1)*dzp
+ DO j=2,nhor-1
+ xcu=x(1)+(j-1)*dxp
+ jcar=1+INT(0.5D0+ &
+ bplina(nx,nz,x,z,rhoplt,xcu,zcu) &
+ *5.D0/density_scale)
+ jcar=MIN(jcar,25)
+ ifun(j)=icar(jcar)
+ END DO
+ IF(MOD(i-1,izsc)==0) THEN
+ ntz=(i+izsc-1)/izsc
+ ifun(nhor)='+'
+ WRITE(*,'(1X,F6.2,1X,A,120A)') zco(ntz),'+',ifun(2:nhor)
+ ELSE
+ ifun(nhor)='i'
+ WRITE(*,'(8X,A,120A)') 'i',ifun(2:nhor)
+ END IF
END DO
IF(MOD(nver-1,izsc)==0) THEN
- WRITE(*,'(1X,F6.2,1X,121A)') zco(ntkz),ibor(1:nhor)
+ WRITE(*,'(1X,F6.2,1X,121A)') zco(ntkz),ibor(1:nhor)
ELSE
- WRITE(*,'(8X,A,120A)') '-',ibor(2:nhor)
+ WRITE(*,'(8X,A,120A)') '-',ibor(2:nhor)
END IF
WRITE(*,'(A,12(F6.2,4X),F6.2)') ' x= ',(xco(i),i=1,ntkx)
END SUBROUTINE plot_density
@@ -172,12 +172,12 @@ CONTAINS
!
ff=0.D0
DO i=1,n-1
- IF(xar(i)<=xcu.AND.xcu<=xar(i+1)) GO TO 10
+ IF(xar(i)<=xcu.AND.xcu<=xar(i+1)) GO TO 10
END DO
RETURN
10 icu=i
DO i=1,m-1
- IF(zar(i)<=zcu.AND.zcu<=zar(i+1)) GO TO 20
+ IF(zar(i)<=zcu.AND.zcu<=zar(i+1)) GO TO 20
END DO
RETURN
20 jcu=i
@@ -201,61 +201,61 @@ CONTAINS
yy=y-cmtot(2)
zz=z-cmtot(3)
DO nst=1,nstmax
- IF(node(nst)/=mpi_myproc) CYCLE
- pst=psi(:,:,:,:,localindex(nst))
- IF(TFFT) THEN
- CALL cdervx(pst,psx)
- CALL cdervy(pst,psy)
+ IF(node(nst)/=mpi_myproc) CYCLE
+ pst=psi(:,:,:,:,localindex(nst))
+ IF(TFFT) THEN
+ CALL cdervx(pst,psx)
+ CALL cdervy(pst,psy)
- CALL cdervz(pst,psz)
- CALL laplace(pst,psw)
- ELSE
- CALL cmulx(der1x,pst,psx,0)
- CALL cmuly(der1y,pst,psy,0)
- CALL cmulz(der1z,pst,psz,0)
- CALL cmulx(der2x,pst,psw,0)
- CALL cmuly(der2y,pst,psw,1)
- CALL cmulz(der2z,pst,psw,1)
- ENDIF
- cc=0.D0
- ss=0.D0
- kin=0.D0
- xpar=0.D0
- DO iz=1,nz
- izz=nz-iz+1
- DO iy=1,ny
- iyy=ny-iy+1
- DO ix=1,nx
- ixx=nx-ix+1
- DO is=1,2
- rp=REAL(pst(ix,iy,iz,is))
- ip=AIMAG(pst(ix,iy,iz,is))
- cc(1)=cc(1)+ &
- rp*(yy(iy)*AIMAG(psz(ix,iy,iz,is))-zz(iz)*AIMAG(psy(ix,iy,iz,is))) &
- +ip*(zz(iz)*REAL(psy(ix,iy,iz,is))-yy(iy)*REAL(psz(ix,iy,iz,is)))
- cc(2)=cc(2)+ &
- rp*(zz(iz)*AIMAG(psx(ix,iy,iz,is))-xx(ix)*AIMAG(psz(ix,iy,iz,is))) &
- +ip*(xx(ix)*REAL(psz(ix,iy,iz,is))-zz(iz)*REAL(psx(ix,iy,iz,is)))
- cc(3)=cc(3)+ &
- rp*(xx(ix)*AIMAG(psy(ix,iy,iz,is))-yy(iy)*AIMAG(psx(ix,iy,iz,is))) &
- +ip*(yy(iy)*REAL(psx(ix,iy,iz,is))-xx(ix)*REAL(psy(ix,iy,iz,is)))
- kin=kin-rp*REAL(psw(ix,iy,iz,is))-ip*AIMAG(psw(ix,iy,iz,is))
- xpar=xpar+REAL(pst(ix,iy,iz,is))*REAL(pst(ixx,iyy,izz,is)) &
- +AIMAG(pst(ix,iy,iz,is))*AIMAG(pst(ixx,iyy,izz,is))
- END DO
- ss(1)=ss(1) + CONJG(pst(ix,iy,iz,1))*pst(ix,iy,iz,2) &
- + CONJG(pst(ix,iy,iz,2))*pst(ix,iy,iz,1)
- ss(2)=ss(2) + CONJG(pst(ix,iy,iz,1))*pst(ix,iy,iz,2)*CMPLX(0.D0,-1.D0,db) &
- + CONJG(pst(ix,iy,iz,2))*pst(ix,iy,iz,1)*CMPLX(0.D0,1.D0,db)
- ss(3)=ss(3) + CONJG(pst(ix,iy,iz,1))*pst(ix,iy,iz,1) &
- - CONJG(pst(ix,iy,iz,2))*pst(ix,iy,iz,2)
+ CALL cdervz(pst,psz)
+ CALL laplace(pst,psw)
+ ELSE
+ CALL cmulx(der1x,pst,psx,0)
+ CALL cmuly(der1y,pst,psy,0)
+ CALL cmulz(der1z,pst,psz,0)
+ CALL cmulx(der2x,pst,psw,0)
+ CALL cmuly(der2y,pst,psw,1)
+ CALL cmulz(der2z,pst,psw,1)
+ ENDIF
+ cc=0.D0
+ ss=0.D0
+ kin=0.D0
+ xpar=0.D0
+ DO iz=1,nz
+ izz=nz-iz+1
+ DO iy=1,ny
+ iyy=ny-iy+1
+ DO ix=1,nx
+ ixx=nx-ix+1
+ DO is=1,2
+ rp=REAL(pst(ix,iy,iz,is))
+ ip=AIMAG(pst(ix,iy,iz,is))
+ cc(1)=cc(1)+ &
+ rp*(yy(iy)*AIMAG(psz(ix,iy,iz,is))-zz(iz)*AIMAG(psy(ix,iy,iz,is))) &
+ +ip*(zz(iz)*REAL(psy(ix,iy,iz,is))-yy(iy)*REAL(psz(ix,iy,iz,is)))
+ cc(2)=cc(2)+ &
+ rp*(zz(iz)*AIMAG(psx(ix,iy,iz,is))-xx(ix)*AIMAG(psz(ix,iy,iz,is))) &
+ +ip*(xx(ix)*REAL(psz(ix,iy,iz,is))-zz(iz)*REAL(psx(ix,iy,iz,is)))
+ cc(3)=cc(3)+ &
+ rp*(xx(ix)*AIMAG(psy(ix,iy,iz,is))-yy(iy)*AIMAG(psx(ix,iy,iz,is))) &
+ +ip*(yy(iy)*REAL(psx(ix,iy,iz,is))-xx(ix)*REAL(psy(ix,iy,iz,is)))
+ kin=kin-rp*REAL(psw(ix,iy,iz,is))-ip*AIMAG(psw(ix,iy,iz,is))
+ xpar=xpar+REAL(pst(ix,iy,iz,is))*REAL(pst(ixx,iyy,izz,is)) &
+ +AIMAG(pst(ix,iy,iz,is))*AIMAG(pst(ixx,iyy,izz,is))
+ END DO
+ ss(1)=ss(1) + CONJG(pst(ix,iy,iz,1))*pst(ix,iy,iz,2) &
+ + CONJG(pst(ix,iy,iz,2))*pst(ix,iy,iz,1)
+ ss(2)=ss(2) + CONJG(pst(ix,iy,iz,1))*pst(ix,iy,iz,2)*CMPLX(0.D0,-1.D0,db) &
+ + CONJG(pst(ix,iy,iz,2))*pst(ix,iy,iz,1)*CMPLX(0.D0,1.D0,db)
+ ss(3)=ss(3) + CONJG(pst(ix,iy,iz,1))*pst(ix,iy,iz,1) &
+ - CONJG(pst(ix,iy,iz,2))*pst(ix,iy,iz,2)
+ ENDDO
ENDDO
- ENDDO
- END DO
- sp_spin(:,nst)=0.5D0*wxyz*ss(:)
- sp_orbital(:,nst)=wxyz*cc
- sp_kinetic(nst)=wxyz*f%h2m(isospin(nst))*kin
- sp_parity(nst)=wxyz*xpar
+ END DO
+ sp_spin(:,nst)=0.5D0*wxyz*ss(:)
+ sp_orbital(:,nst)=wxyz*cc
+ sp_kinetic(nst)=wxyz*f%h2m(isospin(nst))*kin
+ sp_parity(nst)=wxyz*xpar
END DO
END SUBROUTINE sp_properties
!************************************************************
diff --git a/Code/levels.f90 b/Code/levels.f90
index e01ec333d7b58ed9029286a82abd7e94ae3b1bee..eb63a65071b026cf8a378f1707b270a60df356d0 100644
--- a/Code/levels.f90
+++ b/Code/levels.f90
@@ -22,7 +22,7 @@ CONTAINS
isospin(nstmax))
isospin(1:npsi(1))=1
isospin(npmin(2):npsi(2))=2
- END SUBROUTINE alloc_levels
+ END SUBROUTINE alloc_levels
!************************************************************
SUBROUTINE cdervx(psin,d1psout,d2psout)
COMPLEX(db), INTENT(IN) :: psin(:,:,:,:)
diff --git a/Code/main3d.f90 b/Code/main3d.f90
index 7b7cb30ab8242fc5cb58ebf8b18592203b7a4008..bf825ac0acf9c003b5fda163564b1b98f914dc3e 100644
--- a/Code/main3d.f90
+++ b/Code/main3d.f90
@@ -37,40 +37,40 @@ PROGRAM tdhf3d
tdynamic=imode==2
IF(tmpi.AND.tstatic) STOP 'MPI not implemented for static mode'
IF(.NOT.(tstatic.OR.tdynamic)) THEN
- IF(wflag) WRITE(*,*) 'Illegal value for imode:',imode
- STOP
+ IF(wflag) WRITE(*,*) 'Illegal value for imode:',imode
+ STOP
END IF
IF(wflag) THEN
- WRITE(*,*) '***** Main parameter input *****'
- IF(tstatic) WRITE(*,*) 'This is a static calculation'
- IF(tdynamic) WRITE(*,*) 'This is a dynamic calculation'
- WRITE(*,'(3(A16,I5))') 'Print interval:',mprint, &
- 'Plot interval:',mplot,'Save interval:',mrest
- WRITE(*,'(A,F7.4)') 'Radius parameter r0=',r0
- IF(trestart) THEN
- WRITE(*,*) 'The calculation is being restarted from file ',wffile
- ELSE
- SELECT CASE(nof)
+ WRITE(*,*) '***** Main parameter input *****'
+ IF(tstatic) WRITE(*,*) 'This is a static calculation'
+ IF(tdynamic) WRITE(*,*) 'This is a dynamic calculation'
+ WRITE(*,'(3(A16,I5))') 'Print interval:',mprint, &
+ 'Plot interval:',mplot,'Save interval:',mrest
+ WRITE(*,'(A,F7.4)') 'Radius parameter r0=',r0
+ IF(trestart) THEN
+ WRITE(*,*) 'The calculation is being restarted from file ',wffile
+ ELSE
+ SELECT CASE(nof)
CASE(0)
- IF(tdynamic) STOP &
- 'Harmonic oscillator initialization not allowed for dynamic case'
- WRITE(*,*) 'Harmonic oscillator initialization'
+ IF(tdynamic) STOP &
+ 'Harmonic oscillator initialization not allowed for dynamic case'
+ WRITE(*,*) 'Harmonic oscillator initialization'
CASE(1:)
- WRITE(*,'(A,I4,A)') ' Initialization from ',nof,' fragments'
+ WRITE(*,'(A,I4,A)') ' Initialization from ',nof,' fragments'
CASE DEFAULT
- WRITE(*,*) 'User initialization'
- END SELECT
- END IF
- IF(tcoul) THEN
- WRITE(*,*) "Coulomb field is included"
- ELSE
- WRITE(*,*) " Coulomb field *not* included"
- END IF
- WRITE(*,*) "Field output selection: ",writeselect
+ WRITE(*,*) 'User initialization'
+ END SELECT
+ END IF
+ IF(tcoul) THEN
+ WRITE(*,*) "Coulomb field is included"
+ ELSE
+ WRITE(*,*) " Coulomb field *not* included"
+ END IF
+ WRITE(*,*) "Field output selection: ",writeselect
END IF
IF(trestart) THEN
- nofsave=nof ! save whether initial was 2-body
- nof=1 ! restart simulated as fragment input
+ nofsave=nof ! save whether initial was 2-body
+ nof=1 ! restart simulated as fragment input
END IF
IF(.NOT.(tstatic.OR.tdynamic)) STOP 'Illegal value of imode'
!********************************************************************
@@ -84,24 +84,24 @@ PROGRAM tdhf3d
! Step 5: get input for static or dynamic case
!********************************************************************
IF(tstatic) THEN
- CALL getin_static
+ CALL getin_static
ELSE
- CALL getin_dynamic
+ CALL getin_dynamic
END IF
!********************************************************************
! Step 6: get input for static or dynamic case
!********************************************************************
! Determine wave function numbers etc.
IF(nof>mnof) THEN
- STOP 'nof > mnof'
+ STOP 'nof > mnof'
ELSE IF(nof>0) THEN
- CALL getin_fragments
+ CALL getin_fragments
ELSE
- IF(nof==0.AND.tdynamic) STOP 'Dynamic case with nof==0 not allowed'
- npmin(1)=1
- npmin(2)=npsi(1)+1
- npsi(2)=npsi(2)+npsi(1)
- nstmax =npsi(2)
+ IF(nof==0.AND.tdynamic) STOP 'Dynamic case with nof==0 not allowed'
+ npmin(1)=1
+ npmin(2)=npsi(1)+1
+ npsi(2)=npsi(2)+npsi(1)
+ nstmax =npsi(2)
END IF
!********************************************************************
! Step 7: allocate wave functions
@@ -113,15 +113,15 @@ PROGRAM tdhf3d
! Step 8: initialize wave functions
!********************************************************************
IF(nof>0) THEN
- CALL read_fragments
- IF(.NOT.tmpi) THEN
- CALL schmid
- WRITE(*,*) 'Reorthogonalization complete'
- END IF
+ CALL read_fragments
+ IF(.NOT.tmpi) THEN
+ CALL schmid
+ WRITE(*,*) 'Reorthogonalization complete'
+ END IF
ELSEIF(nof==0) THEN
- CALL harmosc
+ CALL harmosc
ELSE
- CALL init_user
+ CALL init_user
END IF
CLOSE(5)
!********************************************************************
@@ -132,16 +132,16 @@ PROGRAM tdhf3d
! Step 10: static or dynamic calculation performed
!********************************************************************
IF(tstatic) THEN
- IF(tmpi .AND. wflag) STOP ' static should not run parallel'
- CALL init_static
- CALL statichf
+ IF(tmpi .AND. wflag) STOP ' static should not run parallel'
+ CALL init_static
+ CALL statichf
ELSE
- !*************************************************************************
- ! Dynamic branch
- !*************************************************************************
- IF(tmpi) CALL mpi_barrier (mpi_comm_world, mpi_ierror)
- IF(trestart) nof=nofsave ! restore 2-body status so analysis is done
- CALL dynamichf
+ !*************************************************************************
+ ! Dynamic branch
+ !*************************************************************************
+ IF(tmpi) CALL mpi_barrier (mpi_comm_world, mpi_ierror)
+ IF(trestart) nof=nofsave ! restore 2-body status so analysis is done
+ CALL dynamichf
ENDIF
CALL finish_mpi
END PROGRAM tdhf3d
diff --git a/Code/meanfield-new.f90 b/Code/meanfield-new.f90
new file mode 100644
index 0000000000000000000000000000000000000000..462b626955f53365d8aa1a21214aa288df96d12d
--- /dev/null
+++ b/Code/meanfield-new.f90
@@ -0,0 +1,247 @@
+Module Meanfield
+ USE Params, ONLY: db,tcoul
+ USE Densities
+ USE Forces
+ USE Grids, ONLY: nx,ny,nz,der1x,der2x,der1y,der2y,der1z,der2z
+ USE Coulomb, ONLY: poisson,wcoul
+ IMPLICIT NONE
+ REAL(db),ALLOCATABLE,DIMENSION(:,:,:,:) :: upot,bmass,divaq
+ REAL(db),ALLOCATABLE,DIMENSION(:,:,:,:,:) :: aq,spot,wlspot,dbmass
+ PRIVATE :: divaq,aq,wlspot,dbmass
+CONTAINS
+ !***********************************************************************
+ SUBROUTINE alloc_fields
+ ALLOCATE(upot(nx,ny,nz,2),bmass(nx,ny,nz,2),divaq(nx,ny,nz,2), &
+ aq(nx,ny,nz,3,2),spot(nx,ny,nz,3,2),wlspot(nx,ny,nz,3,2), &
+ dbmass(nx,ny,nz,3,2))
+ upot=0.D0
+ bmass=0.D0
+ wlspot=0.D0
+ aq=0.D0
+ divaq=0.D0
+ dbmass=0.D0
+ END SUBROUTINE alloc_fields
+ !***********************************************************************
+ SUBROUTINE skyrme
+ USE Trivial, ONLY: rmulx,rmuly,rmulz
+ REAL(db),PARAMETER :: epsilon=1.0d-25
+ REAL(db) :: rotspp,rotspn
+ REAL(db),ALLOCATABLE :: workden(:,:,:,:),workvec(:,:,:,:,:)
+ INTEGER :: ix,iy,iz,ic,iq,icomp
+ ALLOCATE(workden(nx,ny,nz,2),workvec(nx,ny,nz,3,2))
+ ! Step 1: 3-body contribution to upot.
+ DO iq=1,2
+ ic=3-iq
+ upot(:,:,:,iq)=(rho(:,:,:,1)+rho(:,:,:,2))**f%power * &
+ ((b3*(f%power+2.D0)/3.D0-2.D0*b3p/3.D0)*rho(:,:,:,iq) &
+ +b3*(f%power+2.D0)/3.D0*rho(:,:,:,ic) &
+ -(b3p*f%power/3.D0)*(rho(:,:,:,1)**2+rho(:,:,:,2)**2)/ &
+ (rho(:,:,:,1)+rho(:,:,:,2)+epsilon))
+ ENDDO
+ ! Step 2: add divergence of spin-orbit current to upot
+ DO iq=1,2
+ CALL rmulx(der1x,sodens(:,:,:,1,iq),workden(:,:,:,iq),0)
+ CALL rmuly(der1y,sodens(:,:,:,2,iq),workden(:,:,:,iq),1)
+ CALL rmulz(der1z,sodens(:,:,:,3,iq),workden(:,:,:,iq),1)
+ ENDDO
+ DO iq=1,2
+ ic=3-iq
+ upot(:,:,:,iq)=upot(:,:,:,iq) &
+ -(b4+b4p)*workden(:,:,:,iq)-b4*workden(:,:,:,ic)
+ ENDDO
+ ! Step 3: Coulomb potential
+ IF(tcoul) THEN
+ CALL poisson
+ upot(:,:,:,2)=upot(:,:,:,2)+wcoul
+ IF(f%ex/=0) &
+ upot(:,:,:,2)=upot(:,:,:,2)-slate*rho(:,:,:,2)**(1.0D0/3.0D0)
+ ENDIF
+ ! Step 4: remaining terms of upot
+ DO iq=1,2
+ CALL rmulx(der2x,rho(:,:,:,iq),workden(:,:,:,iq),0)
+ CALL rmuly(der2y,rho(:,:,:,iq),workden(:,:,:,iq),1)
+ CALL rmulz(der2z,rho(:,:,:,iq),workden(:,:,:,iq),1)
+ ENDDO
+ DO iq=1,2
+ ic=3-iq
+ upot(:,:,:,iq)=upot(:,:,:,iq)+(b0-b0p)*rho(:,:,:,iq)+b0*rho(:,:,:,ic) &
+ ! t1,t2, and tau-dependent part !
+ +(b1-b1p)*tau(:,:,:,iq)+b1*tau(:,:,:,ic) &
+ ! two-body laplacian*rho-dependent part
+ -(b2-b2p)*workden(:,:,:,iq)-b2*workden(:,:,:,ic)
+ ! Step 5: effective mass
+ bmass(:,:,:,iq)=f%h2m(iq)+(b1-b1p)*rho(:,:,:,iq)+b1*rho(:,:,:,ic)
+ ! Step 6: calculate grad(rho) and wlspot
+ CALL rmulx(der1x,rho(:,:,:,iq),workvec(:,:,:,1,iq),0)
+ CALL rmuly(der1y,rho(:,:,:,iq),workvec(:,:,:,2,iq),0)
+ CALL rmulz(der1z,rho(:,:,:,iq),workvec(:,:,:,3,iq),0)
+ ENDDO
+ DO iq=1,2
+ ic=3-iq
+ wlspot(:,:,:,:,iq)= &
+ (b4+b4p)*workvec(:,:,:,:,iq)+b4*workvec(:,:,:,:,ic)
+ END DO
+ ! Step 7: calculate curl of spin density vector, store in workvec
+ DO iq=1,2
+ CALL rmuly(der1y,sdens(:,:,:,3,iq),workvec(:,:,:,1,iq),0)
+ CALL rmulz(der1z,sdens(:,:,:,2,iq),workvec(:,:,:,1,iq),-1)
+ CALL rmulz(der1z,sdens(:,:,:,1,iq),workvec(:,:,:,2,iq),0)
+ CALL rmulx(der1x,sdens(:,:,:,3,iq),workvec(:,:,:,2,iq),-1)
+ CALL rmulx(der1x,sdens(:,:,:,2,iq),workvec(:,:,:,3,iq),0)
+ CALL rmuly(der1y,sdens(:,:,:,1,iq),workvec(:,:,:,3,iq),-1)
+ ENDDO
+ ! Step 8: calculate A_q vector
+ DO iq=1,2
+ ic=3-iq
+ aq(:,:,:,:,iq)=-2.0D0*(b1-b1p)*current(:,:,:,:,iq) &
+ -2.0D0*b1*current(:,:,:,:,ic) &
+ -(b4+b4p)*workvec(:,:,:,:,iq)-b4*workvec(:,:,:,:,ic)
+ ENDDO
+ ! Step 9: calculate the curl of the current density, stopr in spot
+ DO iq=1,2
+ CALL rmuly(der1y,current(:,:,:,3,iq),spot(:,:,:,1,iq),0)
+ CALL rmulz(der1z,current(:,:,:,2,iq),spot(:,:,:,1,iq),-1)
+ CALL rmulz(der1z,current(:,:,:,1,iq),spot(:,:,:,2,iq),0)
+ CALL rmulx(der1x,current(:,:,:,3,iq),spot(:,:,:,2,iq),-1)
+ CALL rmulx(der1x,current(:,:,:,2,iq),spot(:,:,:,3,iq),0)
+ CALL rmuly(der1y,current(:,:,:,1,iq),spot(:,:,:,3,iq),-1)
+ ENDDO
+ ! Step 10: combine isospin contributions
+ DO icomp=1,3
+ DO iz=1,nz
+ DO iy=1,ny
+ DO ix=1,nx
+ rotspp=spot(ix,iy,iz,icomp,1)
+ rotspn=spot(ix,iy,iz,icomp,2)
+ spot(ix,iy,iz,icomp,1)=-(b4+b4p)*rotspp-b4*rotspn
+ spot(ix,iy,iz,icomp,2)=-(b4+b4p)*rotspn-b4*rotspp
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ! Step 11: calculate divergence of aq in divaq
+ DO iq=1,2
+ CALL rmulx(der1x,aq(:,:,:,1,iq),divaq(:,:,:,iq),0)
+ CALL rmuly(der1y,aq(:,:,:,2,iq),divaq(:,:,:,iq),1)
+ CALL rmulz(der1z,aq(:,:,:,3,iq),divaq(:,:,:,iq),1)
+ ENDDO
+ ! Step 12: calculate the gradient of the effective mass in dbmass
+ DO iq=1,2
+ CALL rmulx(der1x,bmass(:,:,:,iq),dbmass(:,:,:,1,iq),0)
+ CALL rmuly(der1y,bmass(:,:,:,iq),dbmass(:,:,:,2,iq),0)
+ CALL rmulz(der1z,bmass(:,:,:,iq),dbmass(:,:,:,3,iq),0)
+ ENDDO
+ DEALLOCATE(workden,workvec)
+ END SUBROUTINE skyrme
+ !***********************************************************************
+ SUBROUTINE hpsi(iq,eshift,pinn,pout)
+ USE Trivial, ONLY: cmulx, cmuly, cmulz
+ USE Levels, ONLY: cdervx,cdervy,cdervz
+ INTEGER :: iq
+ REAL(db) :: eshift
+ COMPLEX(db),DIMENSION(:,:,:,:) :: pinn,pout
+ INTENT(IN) :: iq,eshift
+ INTENT(INOUT) :: pinn
+ INTENT(OUT) :: pout
+ INTEGER :: is,ic
+ REAL(db) :: sigis
+ COMPLEX(db),ALLOCATABLE,DIMENSION(:,:,:,:) :: pswk,pswk2
+ ALLOCATE(pswk(nx,ny,nz,2),pswk2(nx,ny,nz,2))
+ ! Step 1: non-derivative parts not involving spin
+ DO is=1,2
+ pout(:,:,:,is)=CMPLX(upot(:,:,:,iq)-eshift, &
+ -0D0,db)*pinn(:,:,:,is)
+ ENDDO
+ ! Step 2: the spin-current coupling
+ pout(:,:,:,1)=pout(:,:,:,1) &
+ + CMPLX(spot(:,:,:,1,iq),-spot(:,:,:,2,iq),db) &
+ *pinn(:,:,:,2) + spot(:,:,:,3,iq)*pinn(:,:,:,1)
+ pout(:,:,:,2)=pout(:,:,:,2) &
+ + CMPLX(spot(:,:,:,1,iq),spot(:,:,:,2,iq),db) &
+ *pinn(:,:,:,1) - spot(:,:,:,3,iq)*pinn(:,:,:,2)
+ ! Step 3: derivative terms in x
+ IF(TFFT) THEN
+ CALL cdervx(pinn,pswk,d2psout=pswk2)
+ ELSE
+ CALL cmulx(der1x,pinn,pswk,0)
+ CALL cmulx(der2x,pinn,pswk2,0)
+ ENDIF
+ DO is=1,2
+ ic=3-is
+ sigis=(3-2*is)*0.5D0
+ pout(:,:,:,is)=pout(:,:,:,is) &
+ -CMPLX(dbmass(:,:,:,1,iq),0.5D0*aq(:,:,:,1,iq) &
+ -sigis*wlspot(:,:,:,2,iq),db)*pswk(:,:,:,is) &
+ -sigis*wlspot(:,:,:,3,iq)*pswk(:,:,:,ic) &
+ -bmass(:,:,:,iq)*pswk2(:,:,:,is)
+ ENDDO
+ pswk2(:,:,:,1) = CMPLX(0D0,-0.5D0,db)*&
+ (aq(:,:,:,1,iq)-wlspot(:,:,:,2,iq))*pinn(:,:,:,1)&
+ -0.5D0*wlspot(:,:,:,3,iq)*pinn(:,:,:,2)
+ pswk2(:,:,:,2) = CMPLX(0D0,-0.5D0,db)*&
+ (aq(:,:,:,1,iq)+wlspot(:,:,:,2,iq))*pinn(:,:,:,2)&
+ +0.5D0*wlspot(:,:,:,3,iq)*pinn(:,:,:,1)
+ IF(TFFT) THEN
+ CALL cdervx(pswk2,pswk)
+ ELSE
+ CALL cmulx(der1x,pswk2,pswk,0)
+ ENDIF
+ pout(:,:,:,:)=pout(:,:,:,:) + pswk(:,:,:,:)
+ ! Step 4: derivative terms in y
+ IF(TFFT) THEN
+ CALL cdervy(pinn,pswk,d2psout=pswk2)
+ ELSE
+ CALL cmuly(der1y,pinn,pswk,0)
+ CALL cmuly(der2y,pinn,pswk2,0)
+ ENDIF
+ DO is=1,2
+ ic=3-is
+ sigis=(3-2*is)*0.5D0
+ pout(:,:,:,is)=pout(:,:,:,is) &
+ -CMPLX(dbmass(:,:,:,2,iq),0.5D0*aq(:,:,:,2,iq) &
+ +sigis*wlspot(:,:,:,1,iq),db)*pswk(:,:,:,is) &
+ +CMPLX(0.D0,0.5D0*wlspot(:,:,:,3,iq),db)*pswk(:,:,:,ic) &
+ -bmass(:,:,:,iq)*pswk2(:,:,:,is)
+ ENDDO
+ pswk2(:,:,:,1) = CMPLX(0D0,-0.5D0,db)*&
+ (aq(:,:,:,2,iq)+wlspot(:,:,:,1,iq))*pinn(:,:,:,1)&
+ +CMPLX(0D0,0.5D0,db)*wlspot(:,:,:,3,iq)*pinn(:,:,:,2)
+ pswk2(:,:,:,2) = CMPLX(0D0,-0.5D0,db)*&
+ (aq(:,:,:,2,iq)-wlspot(:,:,:,1,iq))*pinn(:,:,:,2)&
+ +CMPLX(0D0,0.5D0*wlspot(:,:,:,3,iq),db)*pinn(:,:,:,1)
+ IF(TFFT) THEN
+ CALL cdervy(pswk2,pswk)
+ ELSE
+ CALL cmuly(der1y,pswk2,pswk,0)
+ ENDIF
+ pout(:,:,:,:)=pout(:,:,:,:) + pswk(:,:,:,:)
+ ! Step 5: derivative terms in z
+ IF(TFFT) THEN
+ CALL cdervz(pinn,pswk,d2psout=pswk2)
+ ELSE
+ CALL cmulz(der1z,pinn,pswk,0)
+ CALL cmulz(der2z,pinn,pswk2,0)
+ ENDIF
+ DO is=1,2
+ ic=3-is
+ sigis=(3-2*is)*0.5D0
+ pout(:,:,:,is)=pout(:,:,:,is) &
+ -CMPLX(dbmass(:,:,:,3,iq),0.5D0*aq(:,:,:,3,iq),db)*pswk(:,:,:,is) &
+ +CMPLX(sigis*wlspot(:,:,:,1,iq),-0.5D0*wlspot(:,:,:,2,iq),db)* &
+ pswk(:,:,:,ic)-bmass(:,:,:,iq)*pswk2(:,:,:,is)
+ ENDDO
+ pswk2(:,:,:,1) = CMPLX(0D0,-0.5D0,db)*aq(:,:,:,3,iq)*pinn(:,:,:,1)&
+ +CMPLX(0.5D0*wlspot(:,:,:,1,iq),-0.5D0*wlspot(:,:,:,2,iq),db)*pinn(:,:,:,2)
+ pswk2(:,:,:,2) = CMPLX(0D0,-0.5D0,db)*aq(:,:,:,3,iq)*pinn(:,:,:,2)&
+ +CMPLX(-0.5D0*wlspot(:,:,:,1,iq),-0.5D0*wlspot(:,:,:,2,iq),db)*pinn(:,:,:,1)
+ IF(TFFT) THEN
+ CALL cdervz(pswk2,pswk)
+ ELSE
+ CALL cmulz(der1z,pswk2,pswk,0)
+ ENDIF
+ pout(:,:,:,:)=pout(:,:,:,:) + pswk(:,:,:,:)
+ !
+ DEALLOCATE(pswk,pswk2)
+ END SUBROUTINE hpsi
+ !***********************************************************************
+END Module Meanfield
diff --git a/Code/meanfield.f90 b/Code/meanfield.f90
index f489b23b4106bbf3b924cfd2c3aa1a9f01c6ac79..97eb15f3a24143004a6d84920905af3842704229 100644
--- a/Code/meanfield.f90
+++ b/Code/meanfield.f90
@@ -54,7 +54,7 @@ CONTAINS
upot(:,:,:,2)=upot(:,:,:,2)+wcoul
IF(f%ex/=0) &
upot(:,:,:,2)=upot(:,:,:,2)-slate*rho(:,:,:,2)**(1.0D0/3.0D0)
- ENDIF
+ ENDIF
! Step 4: remaining terms of upot
DO iq=1,2
CALL rmulx(der2x,rho(:,:,:,iq),workden(:,:,:,iq),0)
@@ -68,9 +68,9 @@ CONTAINS
+(b1-b1p)*tau(:,:,:,iq)+b1*tau(:,:,:,ic) &
! two-body laplacian*rho-dependent part
-(b2-b2p)*workden(:,:,:,iq)-b2*workden(:,:,:,ic)
- ! Step 5: effective mass
+ ! Step 5: effective mass
bmass(:,:,:,iq)=f%h2m(iq)+(b1-b1p)*rho(:,:,:,iq)+b1*rho(:,:,:,ic)
- ! Step 6: calculate grad(rho) and wlspot
+ ! Step 6: calculate grad(rho) and wlspot
CALL rmulx(der1x,rho(:,:,:,iq),workvec(:,:,:,1,iq),0)
CALL rmuly(der1y,rho(:,:,:,iq),workvec(:,:,:,2,iq),0)
CALL rmulz(der1z,rho(:,:,:,iq),workvec(:,:,:,3,iq),0)
diff --git a/Code/moment.f90 b/Code/moment.f90
index 444d36b974b2f513435f78841ccf039a4035998b..b2eb03b79d760a2a99575eff5cab520d6402fb30 100644
--- a/Code/moment.f90
+++ b/Code/moment.f90
@@ -9,7 +9,7 @@ MODULE Moment
PUBLIC :: pnr,pnrtot,cm,cmtot,pcm,rmstot,beta,gamma, &
moments,moment_print,moment_shortprint
CONTAINS
-!***********************************************************
+ !***********************************************************
SUBROUTINE moments
USE Densities, ONLY: rho,current
INTEGER :: ix,iy,iz,iq
@@ -19,59 +19,59 @@ CONTAINS
cm=0.D0
pcm=0.D0
DO iq=1,2
- DO iz=1,nz
- xx(3)=z(iz)
- DO iy=1,ny
- xx(2)=y(iy)
- DO ix=1,nx
- xx(1)=x(ix)
- pnr(iq)=pnr(iq)+wxyz*rho(ix,iy,iz,iq)
- cm(:,iq)=cm(:,iq)+wxyz*xx*rho(ix,iy,iz,iq)
- IF(tdynamic) pcm(:,iq)=pcm(:,iq)+wxyz*current(ix,iy,iz,:,iq)
+ DO iz=1,nz
+ xx(3)=z(iz)
+ DO iy=1,ny
+ xx(2)=y(iy)
+ DO ix=1,nx
+ xx(1)=x(ix)
+ pnr(iq)=pnr(iq)+wxyz*rho(ix,iy,iz,iq)
+ cm(:,iq)=cm(:,iq)+wxyz*xx*rho(ix,iy,iz,iq)
+ IF(tdynamic) pcm(:,iq)=pcm(:,iq)+wxyz*current(ix,iy,iz,:,iq)
+ ENDDO
ENDDO
- ENDDO
- ENDDO
+ ENDDO
ENDDO
pnrtot=pnr(1)+pnr(2)
cmtot=(cm(:,1)+cm(:,2))/pnrtot
DO iq=1,2
- cm(:,iq)=cm(:,iq)/pnr(iq)
+ cm(:,iq)=cm(:,iq)/pnr(iq)
ENDDO
!***********************************
rms=0.D0
qmat=0.D0
x2m=0.D0
DO iq=1,2
- DO iz=1,nz
- xx(3)=z(iz)-cm(3,iq)
- x2(3)=xx(3)**2
- DO iy=1,ny
- xx(2)=y(iy)-cm(2,iq)
- x2(2)=xx(2)**2
- DO ix=1,nx
- xx(1)=x(ix)-cm(1,iq)
- x2(1)=xx(1)**2
- vol=wxyz*rho(ix,iy,iz,iq)
- rms(iq)=vol*SUM(x2)+rms(iq)
- qmat(1,1,iq)=qmat(1,1,iq)+vol*(x2(1)+x2(1)-x2(2)-x2(3))
- qmat(1,2,iq)=qmat(1,2,iq)+3.D0*vol*xx(1)*xx(2)
- qmat(1,3,iq)=qmat(1,3,iq)+3.D0*vol*xx(1)*xx(3)
- qmat(2,2,iq)=qmat(2,2,iq)+vol*(x2(2)+x2(2)-x2(1)-x2(3))
- qmat(2,3,iq)=qmat(2,3,iq)+3.D0*vol*xx(2)*xx(3)
- qmat(3,3,iq)=qmat(3,3,iq)+vol*(x2(3)+x2(3)-x2(1)-x2(2))
- x2m(:,iq)=vol*x2(:)+x2m(:,iq)
+ DO iz=1,nz
+ xx(3)=z(iz)-cm(3,iq)
+ x2(3)=xx(3)**2
+ DO iy=1,ny
+ xx(2)=y(iy)-cm(2,iq)
+ x2(2)=xx(2)**2
+ DO ix=1,nx
+ xx(1)=x(ix)-cm(1,iq)
+ x2(1)=xx(1)**2
+ vol=wxyz*rho(ix,iy,iz,iq)
+ rms(iq)=vol*SUM(x2)+rms(iq)
+ qmat(1,1,iq)=qmat(1,1,iq)+vol*(x2(1)+x2(1)-x2(2)-x2(3))
+ qmat(1,2,iq)=qmat(1,2,iq)+3.D0*vol*xx(1)*xx(2)
+ qmat(1,3,iq)=qmat(1,3,iq)+3.D0*vol*xx(1)*xx(3)
+ qmat(2,2,iq)=qmat(2,2,iq)+vol*(x2(2)+x2(2)-x2(1)-x2(3))
+ qmat(2,3,iq)=qmat(2,3,iq)+3.D0*vol*xx(2)*xx(3)
+ qmat(3,3,iq)=qmat(3,3,iq)+vol*(x2(3)+x2(3)-x2(1)-x2(2))
+ x2m(:,iq)=vol*x2(:)+x2m(:,iq)
+ ENDDO
ENDDO
ENDDO
- ENDDO
- qmat(2,1,iq)=qmat(1,2,iq)
- qmat(3,1,iq)=qmat(1,3,iq)
- qmat(3,2,iq)=qmat(2,3,iq)
+ qmat(2,1,iq)=qmat(1,2,iq)
+ qmat(3,1,iq)=qmat(1,3,iq)
+ qmat(3,2,iq)=qmat(2,3,iq)
ENDDO
rmstot=SQRT((rms(1)+rms(2))/pnrtot)
rms=SQRT(rms/pnr)
x2mtot=(x2m(:,1)+x2m(:,2))/pnrtot
DO iq=1,2
- x2m(:,iq)=x2m(:,iq)/pnr(iq)
+ x2m(:,iq)=x2m(:,iq)/pnr(iq)
ENDDO
qmtot=qmat(:,:,1)+qmat(:,:,2)
IF(printnow.AND.wflag) WRITE(*,'(/A)') 'Cartesian quadrupole tensor,&
@@ -85,15 +85,15 @@ CONTAINS
beta=SQRT(beta20tot**2+2.0*beta22tot**2)
gamma=ABS(ATAN2(SQRT(2.0)*beta22tot,beta20tot)*180.0D0/PI)
IF(gamma>120.D0) THEN
- gamma=gamma-120.D0
+ gamma=gamma-120.D0
ELSEIF(gamma>60.D0) THEN
- gamma=120.D0-gamma
+ gamma=120.D0-gamma
ENDIF
IF(printnow.AND.wflag) WRITE(*,'(4(A,F8.4)/)') &
' Beta20: ',beta20tot,' Beta22: ',beta22tot,' Beta: ',beta, &
' Gamma: ',gamma
END SUBROUTINE moments
-!***********************************************************
+ !***********************************************************
SUBROUTINE moment_shortprint
OPEN(unit=scratch,file=monopolesfile,POSITION='APPEND')
WRITE(scratch,'(4F10.2,E14.5)') time,rms,rmstot,rms(1)-rms(2)
@@ -112,8 +112,8 @@ CONTAINS
WRITE(*,'(a,2f12.4,1p,4e12.4,3e15.7)') ' Total: ',pnrtot,rmstot, &
q20tot,x2mtot,cmtot
DO iq=1,2
- WRITE(*,'(a,2f12.4,1p,4e12.4,3e15.7)') name(iq),pnr(iq),rms(iq),q20(iq), &
- x2m(:,iq),cm(:,iq)
+ WRITE(*,'(a,2f12.4,1p,4e12.4,3e15.7)') name(iq),pnr(iq),rms(iq),q20(iq), &
+ x2m(:,iq),cm(:,iq)
ENDDO
END SUBROUTINE moment_print
!***********************************************************
@@ -129,8 +129,8 @@ CONTAINS
q_vec=q_mat
IF(info/=0) STOP 'Quadrupole diagonalization failed'
IF(printnow.AND.wflag) THEN
- WRITE(*,'(1X,A,3(F10.2,''('',3F8.4,'')''))') &
- title,(REAL(q_eig(i)),REAL(q_vec(:,i)),i=3,1,-1)
+ WRITE(*,'(1X,A,3(F10.2,''('',3F8.4,'')''))') &
+ title,(REAL(q_eig(i)),REAL(q_vec(:,i)),i=3,1,-1)
ENDIF
q20x=SQRT(5.D0/(16.D0*pi))*REAL(q_eig(3))
q22x=SQRT(5.D0/(96.D0*pi))*(REAL(q_eig(2))-REAL(q_eig(1)))
diff --git a/Code/pairs.f90 b/Code/pairs.f90
index bd10919ecd7dfc66f697f9b64cb03d116c50cf8d..bb897e4732535ea8a51cb01bc05532f87fc90c68 100644
--- a/Code/pairs.f90
+++ b/Code/pairs.f90
@@ -22,20 +22,20 @@ CONTAINS
CALL pairgap
! solve pairing problem for each isospin iq
DO iq=1,2
- IF(iq==2) THEN
- particle_number=charge_number
- ELSE
- particle_number=mass_number-charge_number
- ENDIF
- CALL pairdn(particle_number)
+ IF(iq==2) THEN
+ particle_number=charge_number
+ ELSE
+ particle_number=mass_number-charge_number
+ ENDIF
+ CALL pairdn(particle_number)
ENDDO
! print pairing information
IF(printnow.AND.wflag) THEN
- WRITE(*,'(/7x,a)') ' e_ferm e_pair aver_gap aver_force '
- DO iq=1,2
- WRITE(*,'(a,i2,a,4(1pg12.4))') 'iq=',iq,': ',eferm(iq) , &
- epair(iq) ,avdelt(iq),avg(iq)
- ENDDO
+ WRITE(*,'(/7x,a)') ' e_ferm e_pair aver_gap aver_force '
+ DO iq=1,2
+ WRITE(*,'(a,i2,a,4(1pg12.4))') 'iq=',iq,': ',eferm(iq) , &
+ epair(iq) ,avdelt(iq),avg(iq)
+ ENDDO
ENDIF
END SUBROUTINE pair
!***********************************************************************
@@ -51,40 +51,40 @@ CONTAINS
REAL(db) :: work(nx,ny,nz)
! constant gap in early stages of iteration
IF(iter<=itrsin) THEN
- deltaf=11.2/SQRT(mass_number)
- RETURN
+ deltaf=11.2/SQRT(mass_number)
+ RETURN
ENDIF
! now the detailed gaps:
DO iqq=1,2
- work=0.D0
- ! accumulate new pair-density
- DO nst=npmin(iqq),npsi(iqq)
- DO is=1,2
- work=work+SQRT(MAX(wocc(nst)-wocc(nst)**2,smallp))*0.5* &
- (REAL(psi(:,:,:,is,nst))**2+AIMAG(psi(:,:,:,is,nst))**2)
- ENDDO
- ENDDO
- ! determine pairing strength
- IF(iqq==2) THEN
- v0act=p%v0prot
- ELSE
- v0act=p%v0neut
- ENDIF
- ! now multiply with strength to obtain local pair-potential
- IF(ipair==6) THEN
- work=v0act*work*(1D0-(rho(:,:,:,1)+rho(:,:,:,2))/p%rho0pr)
- ELSE
- work=v0act*work
- END IF
- ! finally compute the actual gaps as s.p. expectation values with
- ! the pair potential
- DO nst=npmin(iqq),npsi(iqq)
- deltaf(nst)=0.0D0
- DO is=1,2
- deltaf(nst)=deltaf(nst)+wxyz*SUM(work* &
- (REAL(psi(:,:,:,is,nst))**2+AIMAG(psi(:,:,:,is,nst))**2))
- ENDDO
- ENDDO
+ work=0.D0
+ ! accumulate new pair-density
+ DO nst=npmin(iqq),npsi(iqq)
+ DO is=1,2
+ work=work+SQRT(MAX(wocc(nst)-wocc(nst)**2,smallp))*0.5* &
+ (REAL(psi(:,:,:,is,nst))**2+AIMAG(psi(:,:,:,is,nst))**2)
+ ENDDO
+ ENDDO
+ ! determine pairing strength
+ IF(iqq==2) THEN
+ v0act=p%v0prot
+ ELSE
+ v0act=p%v0neut
+ ENDIF
+ ! now multiply with strength to obtain local pair-potential
+ IF(ipair==6) THEN
+ work=v0act*work*(1D0-(rho(:,:,:,1)+rho(:,:,:,2))/p%rho0pr)
+ ELSE
+ work=v0act*work
+ END IF
+ ! finally compute the actual gaps as s.p. expectation values with
+ ! the pair potential
+ DO nst=npmin(iqq),npsi(iqq)
+ deltaf(nst)=0.0D0
+ DO is=1,2
+ deltaf(nst)=deltaf(nst)+wxyz*SUM(work* &
+ (REAL(psi(:,:,:,is,nst))**2+AIMAG(psi(:,:,:,is,nst))**2))
+ ENDDO
+ ENDDO
ENDDO
END SUBROUTINE pairgap
!***********************************************************************
@@ -108,13 +108,13 @@ CONTAINS
sumduv=0.0D0
sumepa=0.0D0
DO na=npmin(iq),npsi(iq)
- edif=sp_energy(na)-eferm(iq)
- equasi=SQRT(edif*edif+deltaf(na)**2)
- v2=0.5D0-0.5D0*edif/equasi
- vol=0.5D0*SQRT(MAX(v2-v2*v2,xsmall))
- sumuv=vol+sumuv
- sumduv=vol*deltaf(na)+sumduv
- sumepa=0.5D0*deltaf(na)**2/equasi+sumepa
+ edif=sp_energy(na)-eferm(iq)
+ equasi=SQRT(edif*edif+deltaf(na)**2)
+ v2=0.5D0-0.5D0*edif/equasi
+ vol=0.5D0*SQRT(MAX(v2-v2*v2,xsmall))
+ sumuv=vol+sumuv
+ sumduv=vol*deltaf(na)+sumduv
+ sumepa=0.5D0*deltaf(na)**2/equasi+sumepa
ENDDO
sumuv=MAX(sumuv,xsmall)
avdelt(iq)=sumduv/sumuv
@@ -145,83 +145,83 @@ CONTAINS
rd=rb-ra
re=rd
IF((rfa>0.D0) .AND.(rfb>0.D0)) THEN
- IF(rfa>rfb) THEN
- rb=1.0
- ELSE
- rb=0.D0
- ENDIF
+ IF(rfa>rfb) THEN
+ rb=1.0
+ ELSE
+ rb=0.D0
+ ENDIF
ELSEIF((rfa<0.D0) .AND.(rfb<0.D0)) THEN
- IF(rfa>rfb) THEN
- rb=0.D0
- ELSE
- rb=1.0
- ENDIF
+ IF(rfa>rfb) THEN
+ rb=0.D0
+ ELSE
+ rb=1.0
+ ENDIF
ELSE
- iteration: DO ii=1,iitmax
- IF(((rfb>0.D0) .AND.(rfc>0.D0)).OR.((rfb<0.D0) &
- .AND.(rfc<0.D0))) THEN
- rc=ra
- rfc=rfa
- rd=rb-ra
- re=rd
- ENDIF
- IF(ABS(rfc)<ABS(rfb)) THEN
- ra=rb
- rb=rc
- rc=ra
- rfa=rfb
- rfb=rfc
- rfc=rfa
- ENDIF
- ! convergence check
- rtol1=2.0*reps*ABS(rb)+0.5*rtol
- rxm=0.5 *(rc-rb)
- IF((ABS(rxm) <=rtol1).OR.((ABS(rfb)==0.D0))) THEN
- res=rb
- CALL bcs_occupation(res,rfa)
- RETURN
- ENDIF
- IF((ABS(re)>=rtol1).OR.(ABS(rfa)>ABS(rfb))) &
- THEN
- rs=rfb/rfa
- IF(ra==rc) THEN
- rp=2.0*rxm*rs
- rq=1.0-rs
- ELSE
- rq=rfa/rfc
- rr=rfb/rfc
- rp=rs *(2.0*rxm*rq *(rq-rr)-(rb-ra) *(rr- &
- 1.0))
- rq=(rq-1.0) *(rr-1.0) *(rs-1.0)
+ iteration: DO ii=1,iitmax
+ IF(((rfb>0.D0) .AND.(rfc>0.D0)).OR.((rfb<0.D0) &
+ .AND.(rfc<0.D0))) THEN
+ rc=ra
+ rfc=rfa
+ rd=rb-ra
+ re=rd
+ ENDIF
+ IF(ABS(rfc)<ABS(rfb)) THEN
+ ra=rb
+ rb=rc
+ rc=ra
+ rfa=rfb
+ rfb=rfc
+ rfc=rfa
ENDIF
- IF(rp>0.D0) THEN
- rq=-rq
+ ! convergence check
+ rtol1=2.0*reps*ABS(rb)+0.5*rtol
+ rxm=0.5 *(rc-rb)
+ IF((ABS(rxm) <=rtol1).OR.((ABS(rfb)==0.D0))) THEN
+ res=rb
+ CALL bcs_occupation(res,rfa)
+ RETURN
ENDIF
- rp=ABS(rp)
- rmin1=3.0*rxm*rq-ABS(rtol1*rq)
- rmin2=ABS(rq*re)
- IF(2.0*rp<MIN(rmin1,rmin2)) THEN
- re=rd
- rd=rp/rq
+ IF((ABS(re)>=rtol1).OR.(ABS(rfa)>ABS(rfb))) &
+ THEN
+ rs=rfb/rfa
+ IF(ra==rc) THEN
+ rp=2.0*rxm*rs
+ rq=1.0-rs
+ ELSE
+ rq=rfa/rfc
+ rr=rfb/rfc
+ rp=rs *(2.0*rxm*rq *(rq-rr)-(rb-ra) *(rr- &
+ 1.0))
+ rq=(rq-1.0) *(rr-1.0) *(rs-1.0)
+ ENDIF
+ IF(rp>0.D0) THEN
+ rq=-rq
+ ENDIF
+ rp=ABS(rp)
+ rmin1=3.0*rxm*rq-ABS(rtol1*rq)
+ rmin2=ABS(rq*re)
+ IF(2.0*rp<MIN(rmin1,rmin2)) THEN
+ re=rd
+ rd=rp/rq
+ ELSE
+ rd=rxm
+ re=rd
+ ENDIF
+ ELSE
+ rd=rxm
+ re=rd
+ ENDIF
+ ra=rb
+ rfa=rfb
+ IF(ABS(rd) >rtol1) THEN
+ rb=rb+rd
ELSE
- rd=rxm
- re=rd
+ rb=rb+SIGN(rtol1,rxm)
ENDIF
- ELSE
- rd=rxm
- re=rd
- ENDIF
- ra=rb
- rfa=rfb
- IF(ABS(rd) >rtol1) THEN
- rb=rb+rd
- ELSE
- rb=rb+SIGN(rtol1,rxm)
- ENDIF
- CALL bcs_occupation(rb,rfb)
- rfb=particle_number-rfb
- ENDDO iteration
- STOP 'No solution found in pairing iterations'
+ CALL bcs_occupation(rb,rfb)
+ rfb=particle_number-rfb
+ ENDDO iteration
+ STOP 'No solution found in pairing iterations'
ENDIF
END FUNCTION rbrent
!***********************************************************************
@@ -235,10 +235,10 @@ CONTAINS
REAL(db) :: edif
bcs_partnum=0.0
DO k=npmin(iq),npsi(iq)
- edif=sp_energy(k)-efermi
- wocc(k)=0.5D0 *(1.0D0-edif/SQRT(edif**2+deltaf(k)**2))
- wocc(k)=MIN(MAX(wocc(k),smal),1.D0-smal)
- bcs_partnum=bcs_partnum+wocc(k)
+ edif=sp_energy(k)-efermi
+ wocc(k)=0.5D0 *(1.0D0-edif/SQRT(edif**2+deltaf(k)**2))
+ wocc(k)=MIN(MAX(wocc(k),smal),1.D0-smal)
+ bcs_partnum=bcs_partnum+wocc(k)
ENDDO
END SUBROUTINE bcs_occupation
END MODULE Pairs
diff --git a/Code/parallel.f90 b/Code/parallel.f90
index 89bef772165a8c65cfe42be77930a7877a853ffb..9038219d659b3a39195ddcb3a93198e51be15d72 100644
--- a/Code/parallel.f90
+++ b/Code/parallel.f90
@@ -9,11 +9,11 @@ MODULE Parallel
INTEGER, ALLOCATABLE :: node(:),localindex(:),globalindex(:)
INTEGER :: mpi_nprocs,mpi_ierror,mpi_myproc
CONTAINS
-!***********************************************************************
+ !***********************************************************************
SUBROUTINE alloc_nodes
ALLOCATE(node(nstmax),localindex(nstmax),globalindex(nstmax))
END SUBROUTINE alloc_nodes
-!***********************************************************************
+ !***********************************************************************
SUBROUTINE init_all_mpi
CALL mpi_init(mpi_ierror)
CALL mpi_comm_size(mpi_comm_world,mpi_nprocs,mpi_ierror)
@@ -22,7 +22,7 @@ CONTAINS
IF(wflag) WRITE(*,'(a,i5)') ' number of nodes=',mpi_nprocs
CALL mpi_barrier (mpi_comm_world, mpi_ierror)
END SUBROUTINE init_all_mpi
-!***********************************************************************
+ !***********************************************************************
SUBROUTINE associate_nodes
INTEGER :: ncount,nst,ip,iloc
ncount=0
@@ -39,13 +39,13 @@ CONTAINS
ENDIF
ENDDO
DO ip=0,mpi_nprocs-1
- iloc=0
- DO nst=1,nstmax
- IF(node(nst)==ip) THEN
- iloc=iloc+1
- localindex(nst)=iloc
- ENDIF
- ENDDO
+ iloc=0
+ DO nst=1,nstmax
+ IF(node(nst)==ip) THEN
+ iloc=iloc+1
+ localindex(nst)=iloc
+ ENDIF
+ ENDDO
ENDDO
IF(wflag) THEN
WRITE(*,'(A/(1X,20I4))') &
@@ -53,7 +53,7 @@ CONTAINS
ENDIF
CALL mpi_barrier (mpi_comm_world, mpi_ierror)
END SUBROUTINE associate_nodes
-!***********************************************************************
+ !***********************************************************************
SUBROUTINE collect_densities
USE Densities, ONLY : rho,tau,current,sodens,sdens
REAL(db) :: tmp_rho(nx,ny,nz,2),tmp_current(nx,ny,nz,3,2)
@@ -78,7 +78,7 @@ CONTAINS
CALL mpi_barrier (mpi_comm_world,mpi_ierror)
RETURN
END SUBROUTINE collect_densities
-!***********************************************************************
+ !***********************************************************************
SUBROUTINE collect_sp_properties
REAL(db) :: tmpgat(nstmax),tmpgat3(3,nstmax)
CALL mpi_barrier (mpi_comm_world,mpi_ierror)
@@ -108,7 +108,7 @@ CONTAINS
sp_parity=tmpgat
CALL mpi_barrier (mpi_comm_world,mpi_ierror)
END SUBROUTINE collect_sp_properties
-!***********************************************************************
+ !***********************************************************************
SUBROUTINE finish_mpi
INTEGER :: ierr
CALL mpi_finalize(ierr)
diff --git a/Code/static.f90 b/Code/static.f90
index 60f085ff1b25be0a119fffd7bde7e973ad321a79..f5f7230f414f2c7c5fe14dcef610451948e95ee5 100644
--- a/Code/static.f90
+++ b/Code/static.f90
@@ -23,47 +23,47 @@ CONTAINS
npsi=0
READ(5,static)
IF(nof<=0) THEN
- IF(npsi(1)==0) THEN
- IF(ipair==0.OR.nof<0) THEN
- npsi(1)=nneut
- ELSE
- npsi(1)=NINT(nneut+1.65*FLOAT(nneut)**0.666667D0)
- ENDIF
- ENDIF
- IF(npsi(2)==0) THEN
- IF(ipair==0.OR.nof<0) THEN
- npsi(2)=nprot
- ELSE
- npsi(2)=NINT(nprot+1.65*FLOAT(nprot)**0.666667D0)
- ENDIF
- ENDIF
- IF(nneut>npsi(1).OR.nprot>npsi(2)) &
- STOP 'Particle & state numbers in conflict'
- nstmax=npsi(1)+npsi(2)
- charge_number=nprot
- mass_number=nneut+nprot
+ IF(npsi(1)==0) THEN
+ IF(ipair==0.OR.nof<0) THEN
+ npsi(1)=nneut
+ ELSE
+ npsi(1)=NINT(nneut+1.65*FLOAT(nneut)**0.666667D0)
+ ENDIF
+ ENDIF
+ IF(npsi(2)==0) THEN
+ IF(ipair==0.OR.nof<0) THEN
+ npsi(2)=nprot
+ ELSE
+ npsi(2)=NINT(nprot+1.65*FLOAT(nprot)**0.666667D0)
+ ENDIF
+ ENDIF
+ IF(nneut>npsi(1).OR.nprot>npsi(2)) &
+ STOP 'Particle & state numbers in conflict'
+ nstmax=npsi(1)+npsi(2)
+ charge_number=nprot
+ mass_number=nneut+nprot
END IF
END SUBROUTINE getin_static
!*************************************************************************
SUBROUTINE init_static
IF(wflag) THEN
- WRITE(*,*) '***** Parameters for static calculation *****'
- WRITE(*,"(3(A,I4))") ' Neutrons:',nneut,' in ',npsi(1),' levels'
- WRITE(*,"(3(A,I4))") ' Protons :',nprot,' in ',npsi(2)-npsi(1),' levels'
- WRITE(*,"(A,I6)") " Maximum number of iterations: ",maxiter
- WRITE(*,"(2(A,G12.5))") ' Damping coefficient:',x0dmp, &
- " Damping energy scale: ",e0dmp
- WRITE(*,"(A,1PE12.4)") " Convergence limit: ",serr
- ! initialize *.res files
- CALL start_protocol(converfile, &
- '# Iter Energy d_Energy h**2 h*h rms &
- &beta2 gamma')
- CALL start_protocol(dipolesfile, &
- '# Iter c.m. x-y-z Isovector&
- &dipoles x-y-z')
- CALL start_protocol(spinfile, &
- '# Iter Lx Ly Lz Sx Sy &
- &Sz Jx Jy Jz')
+ WRITE(*,*) '***** Parameters for static calculation *****'
+ WRITE(*,"(3(A,I4))") ' Neutrons:',nneut,' in ',npsi(1),' levels'
+ WRITE(*,"(3(A,I4))") ' Protons :',nprot,' in ',npsi(2)-npsi(1),' levels'
+ WRITE(*,"(A,I6)") " Maximum number of iterations: ",maxiter
+ WRITE(*,"(2(A,G12.5))") ' Damping coefficient:',x0dmp, &
+ " Damping energy scale: ",e0dmp
+ WRITE(*,"(A,1PE12.4)") " Convergence limit: ",serr
+ ! initialize *.res files
+ CALL start_protocol(converfile, &
+ '# Iter Energy d_Energy h**2 h*h rms &
+ &beta2 gamma')
+ CALL start_protocol(dipolesfile, &
+ '# Iter c.m. x-y-z Isovector&
+ &dipoles x-y-z')
+ CALL start_protocol(spinfile, &
+ '# Iter Lx Ly Lz Sx Sy &
+ &Sz Jx Jy Jz')
ENDIF
! calculate damping matrices
IF(e0dmp>0.0D0) CALL setup_damping(e0dmp)
@@ -82,16 +82,16 @@ CONTAINS
! Step 1: initialization
IF(tdiag) ALLOCATE(hmatr(nstmax,nstmax))
IF(trestart) THEN
- firstiter=iter+1
+ firstiter=iter+1
ELSE
- iter=0
- firstiter=1
- sp_energy=0.0D0
- sp_efluct1=0.0D0
- sp_efluct2=0.D0
- sp_norm=0.0D0
- sumflu=0.D0
- CALL schmid
+ iter=0
+ firstiter=1
+ sp_energy=0.0D0
+ sp_efluct1=0.0D0
+ sp_efluct2=0.D0
+ sp_norm=0.0D0
+ sumflu=0.D0
+ CALL schmid
END IF
! Step 2: calculate densities and mean field
rho=0.0D0
@@ -100,22 +100,22 @@ CONTAINS
sdens=0.0D0
sodens=0.0D0
DO nst=1,nstmax
- CALL add_density(isospin(nst),wocc(nst),psi(:,:,:,:,nst), &
- rho,tau,current,sdens,sodens)
+ CALL add_density(isospin(nst),wocc(nst),psi(:,:,:,:,nst), &
+ rho,tau,current,sdens,sodens)
ENDDO
CALL skyrme
! Step 3: initial gradient step
delesum=0.0D0
sumflu=0.0D0
-!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,denerg) &
-!$OMP SCHEDULE(STATIC) REDUCTION(+: sumflu , delesum)
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,denerg) &
+ !$OMP SCHEDULE(STATIC) REDUCTION(+: sumflu , delesum)
DO nst=1,nstmax
- CALL grstep(nst,isospin(nst),sp_energy(nst),denerg, &
- psi(:,:,:,:,nst))
- sumflu=sumflu+wocc(nst)*sp_efluct1(nst)
- delesum=delesum+wocc(nst)*denerg
+ CALL grstep(nst,isospin(nst),sp_energy(nst),denerg, &
+ psi(:,:,:,:,nst))
+ sumflu=sumflu+wocc(nst)*sp_efluct1(nst)
+ delesum=delesum+wocc(nst)*denerg
ENDDO
-!$OMP END PARALLEL DO
+ !$OMP END PARALLEL DO
! pairing and orthogonalization
IF(ipair/=0) CALL pair
CALL schmid
@@ -124,65 +124,65 @@ CONTAINS
CALL sinfo
! step 4: start static iteration loop
Iteration: DO iter=firstiter,maxiter
- WRITE(*,'(a,i6)') ' Static Iteration No.',iter
- ! Step 5: gradient step
- delesum=0.0D0
- sumflu=0.0D0
-!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,denerg) &
-!$OMP SCHEDULE(STATIC) REDUCTION(+: sumflu , delesum)
- DO nst=1,nstmax
- CALL grstep(nst,isospin(nst),sp_energy(nst),denerg, &
- psi(:,:,:,:,nst))
- sumflu=sumflu+wocc(nst)*sp_efluct1(nst)
- delesum=delesum+wocc(nst)*denerg
- ENDDO
-!$OMP END PARALLEL DO
- ! Step 6: diagonalize if desired
- IF(tdiag.AND.iter>20) THEN
-!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(iq) SCHEDULE(STATIC)
- DO iq=1,2
- CALL diagstep(iq,npsi(iq)-npmin(iq)+1)
- ENDDO
-!$OMP END PARALLEL DO
- ENDIF
- ! Step 7: do pairing and orthogonalization
- IF(ipair/=0) CALL pair
- CALL schmid
- ! Step 8: get new densities and fields with relaxation
- IF(taddnew) THEN
- upot=rho
- bmass=tau
- ENDIF
- rho=0.0D0
- tau=0.0D0
- current=0.0D0
- sdens=0.0D0
- sodens=0.0D0
-!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) &
-!$OMP REDUCTION(+:rho, tau, current, sdens, sodens)
- DO nst=1,nstmax
- CALL add_density(isospin(nst),wocc(nst),psi(:,:,:,:,nst), &
- rho,tau,current,sdens,sodens)
- ENDDO
-!$OMP END PARALLEL DO
- IF(taddnew) THEN
- rho=addnew*rho+addco*upot
- tau=addnew*tau+addco*bmass
- ENDIF
- CALL skyrme
- ! calculate and print information
- IF(mprint>0.AND.MOD(iter,mprint)==0) THEN
- CALL sp_properties
- CALL sinfo
- ENDIF
- ! Step 9: check for convergence, saving wave functions
- IF(sumflu/nstmax<serr.AND.iter>1) THEN
- CALL write_wavefunctions
- EXIT Iteration
- END IF
- IF(MOD(iter,mrest)==0) THEN
- CALL write_wavefunctions
- ENDIF
+ WRITE(*,'(a,i6)') ' Static Iteration No.',iter
+ ! Step 5: gradient step
+ delesum=0.0D0
+ sumflu=0.0D0
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,denerg) &
+ !$OMP SCHEDULE(STATIC) REDUCTION(+: sumflu , delesum)
+ DO nst=1,nstmax
+ CALL grstep(nst,isospin(nst),sp_energy(nst),denerg, &
+ psi(:,:,:,:,nst))
+ sumflu=sumflu+wocc(nst)*sp_efluct1(nst)
+ delesum=delesum+wocc(nst)*denerg
+ ENDDO
+ !$OMP END PARALLEL DO
+ ! Step 6: diagonalize if desired
+ IF(tdiag.AND.iter>20) THEN
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(iq) SCHEDULE(STATIC)
+ DO iq=1,2
+ CALL diagstep(iq,npsi(iq)-npmin(iq)+1)
+ ENDDO
+ !$OMP END PARALLEL DO
+ ENDIF
+ ! Step 7: do pairing and orthogonalization
+ IF(ipair/=0) CALL pair
+ CALL schmid
+ ! Step 8: get new densities and fields with relaxation
+ IF(taddnew) THEN
+ upot=rho
+ bmass=tau
+ ENDIF
+ rho=0.0D0
+ tau=0.0D0
+ current=0.0D0
+ sdens=0.0D0
+ sodens=0.0D0
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) &
+ !$OMP REDUCTION(+:rho, tau, current, sdens, sodens)
+ DO nst=1,nstmax
+ CALL add_density(isospin(nst),wocc(nst),psi(:,:,:,:,nst), &
+ rho,tau,current,sdens,sodens)
+ ENDDO
+ !$OMP END PARALLEL DO
+ IF(taddnew) THEN
+ rho=addnew*rho+addco*upot
+ tau=addnew*tau+addco*bmass
+ ENDIF
+ CALL skyrme
+ ! calculate and print information
+ IF(mprint>0.AND.MOD(iter,mprint)==0) THEN
+ CALL sp_properties
+ CALL sinfo
+ ENDIF
+ ! Step 9: check for convergence, saving wave functions
+ IF(sumflu/nstmax<serr.AND.iter>1) THEN
+ CALL write_wavefunctions
+ EXIT Iteration
+ END IF
+ IF(MOD(iter,mrest)==0) THEN
+ CALL write_wavefunctions
+ ENDIF
END DO Iteration
IF(tdiag) DEALLOCATE(hmatr)
END SUBROUTINE statichf
@@ -225,19 +225,19 @@ CONTAINS
ENDIF
! Step 4: the damping step
IF(e0dmp>0.0D0) THEN
- ps1=ps1 - xnormb*psin
- x0act=x0dmp
- IF(TFFT) THEN
- CALL laplace(ps1,ps2,e0inv=e0dmp)
- ELSE
- CALL cmulz(cdmpz,ps1,ps2,0)
- CALL cmuly(cdmpy,ps2,ps1,0)
- CALL cmulx(cdmpx,ps1,ps2,0)
- x0act=x0act/e0dmp
- ENDIF
- psin=psin - x0act*ps2
+ ps1=ps1 - xnormb*psin
+ x0act=x0dmp
+ IF(TFFT) THEN
+ CALL laplace(ps1,ps2,e0inv=e0dmp)
+ ELSE
+ CALL cmulz(cdmpz,ps1,ps2,0)
+ CALL cmuly(cdmpy,ps2,ps1,0)
+ CALL cmulx(cdmpx,ps1,ps2,0)
+ x0act=x0act/e0dmp
+ ENDIF
+ psin=psin - x0act*ps2
ELSE
- psin=(1.0+x0dmp*xnormb)*psin-x0dmp*ps1
+ psin=(1.0+x0dmp*xnormb)*psin-x0dmp*ps1
ENDIF
! Step 5: energy convergence criterion
enrold=spe
@@ -263,60 +263,60 @@ CONTAINS
REAL(db) :: rwork(2*nlin*nlin+5*nlin+1)
INTEGER :: iwork(5*nlin+3)
INTERFACE
- SUBROUTINE zhbevd( jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, &
- lwork, rwork, lrwork, iwork, liwork, info )
- USE Params, ONLY: db
- CHARACTER(1) :: jobz, uplo
- INTEGER :: info, kd, ldab, ldz, liwork, lrwork, lwork, n, iwork(*)
- DOUBLE PRECISION :: rwork( * ), w( * )
- COMPLEX(8) :: ab( ldab, * ), work( * ), z( ldz, * )
- INTENT(IN) :: jobz,uplo,n,kd,ldab,ldz,lwork,lrwork,liwork
- INTENT(INOUT) :: ab
- INTENT(OUT) :: w,z,work,rwork,iwork,info
- END SUBROUTINE zhbevd
+ SUBROUTINE zhbevd( jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, &
+ lwork, rwork, lrwork, iwork, liwork, info )
+ USE Params, ONLY: db
+ CHARACTER(1) :: jobz, uplo
+ INTEGER :: info, kd, ldab, ldz, liwork, lrwork, lwork, n, iwork(*)
+ DOUBLE PRECISION :: rwork( * ), w( * )
+ COMPLEX(8) :: ab( ldab, * ), work( * ), z( ldz, * )
+ INTENT(IN) :: jobz,uplo,n,kd,ldab,ldz,lwork,lrwork,liwork
+ INTENT(INOUT) :: ab
+ INTENT(OUT) :: w,z,work,rwork,iwork,info
+ END SUBROUTINE zhbevd
END INTERFACE
! Step 1: copy matrix into symmetric storage mode, then diagonalize
noffset=npmin(iq)-1
DO nst=npmin(iq),npsi(iq)
- DO nst2=nst,npsi(iq)
- hmatr_lin(1+nst2-nst,nst-noffset)=&
- 0.5D0*(CONJG(hmatr(nst,nst2))+hmatr(nst2,nst))
- ENDDO
+ DO nst2=nst,npsi(iq)
+ hmatr_lin(1+nst2-nst,nst-noffset)=&
+ 0.5D0*(CONJG(hmatr(nst,nst2))+hmatr(nst2,nst))
+ ENDDO
ENDDO
CALL ZHBEVD('V','L',nlin,nlin,hmatr_lin,nlin+1,eigen,unitary,nstmax, &
cwork,nlin*nlin*2,rwork,2*nlin*nlin+5*nlin+1, &
iwork,5*nlin+3,infoconv)
! Step 2: transform states, replace original ones
IF(tlarge) THEN
- OPEN(scratch,status='scratch',form='unformatted')
- ALLOCATE(psiw(nx,ny,nz,2))
- noffset=npmin(iq)-1
- DO nst=npmin(iq),npsi(iq)
- psiw=CMPLX(0.0D0,0.0D0,db)
- DO nst2=npmin(iq),npsi(iq)
- psiw(:,:,:,:)=psiw(:,:,:,:) &
- + unitary(nst2-noffset,nst-noffset)*psi(:,:,:,:,nst2)
- ENDDO
- WRITE(scratch) psiw
- ENDDO
- REWIND(scratch)
- DO nst=npmin(iq),npsi(iq)
- READ(scratch) psi(:,:,:,:,nst)
- ENDDO
- DEALLOCATE(psiw)
- CLOSE(scratch)
+ OPEN(scratch,status='scratch',form='unformatted')
+ ALLOCATE(psiw(nx,ny,nz,2))
+ noffset=npmin(iq)-1
+ DO nst=npmin(iq),npsi(iq)
+ psiw=CMPLX(0.0D0,0.0D0,db)
+ DO nst2=npmin(iq),npsi(iq)
+ psiw(:,:,:,:)=psiw(:,:,:,:) &
+ + unitary(nst2-noffset,nst-noffset)*psi(:,:,:,:,nst2)
+ ENDDO
+ WRITE(scratch) psiw
+ ENDDO
+ REWIND(scratch)
+ DO nst=npmin(iq),npsi(iq)
+ READ(scratch) psi(:,:,:,:,nst)
+ ENDDO
+ DEALLOCATE(psiw)
+ CLOSE(scratch)
ELSE
- ALLOCATE(ps1(nx,ny,nz,2,nlin))
- noffset=npmin(iq)-1
- ps1=0.0D0
- DO nst=npmin(iq),npsi(iq)
- DO nst2=npmin(iq),npsi(iq)
- ps1(:,:,:,:,nst-noffset)=ps1(:,:,:,:,nst-noffset) &
- + unitary(nst2-noffset,nst-noffset)*psi(:,:,:,:,nst2)
- ENDDO
- ENDDO
- psi(:,:,:,:,npmin(iq):npsi(iq))=ps1
- DEALLOCATE(ps1)
+ ALLOCATE(ps1(nx,ny,nz,2,nlin))
+ noffset=npmin(iq)-1
+ ps1=0.0D0
+ DO nst=npmin(iq),npsi(iq)
+ DO nst2=npmin(iq),npsi(iq)
+ ps1(:,:,:,:,nst-noffset)=ps1(:,:,:,:,nst-noffset) &
+ + unitary(nst2-noffset,nst-noffset)*psi(:,:,:,:,nst2)
+ ENDDO
+ ENDDO
+ psi(:,:,:,:,npmin(iq):npsi(iq))=ps1
+ DEALLOCATE(ps1)
ENDIF
END SUBROUTINE diagstep
!*************************************************************************
@@ -356,22 +356,22 @@ CONTAINS
epair(1),' protons: ',epair(2)
! output densities
IF(mplot/=0) THEN
- IF(MOD(iter,mplot)==0) THEN
- CALL plot_density
- CALL write_densities
- ENDIF
+ IF(MOD(iter,mplot)==0) THEN
+ CALL plot_density
+ CALL write_densities
+ ENDIF
ENDIF
IF(.NOT.wflag) RETURN
! print details of s.p. levels
WRITE(*,'(A)') ' Neutron Single Particle States:',header
DO il=1,nstmax
- IF(il==npmin(2)) THEN
- WRITE(*,'(A)') ' Proton Single Particle States:',header
- END IF
- WRITE(*,'(1X,I3,F4.0,F8.5,2F9.5,F9.6,F8.3,F10.3,3F8.3,3F7.3)') &
- il,sp_parity(il),wocc(il),sp_efluct1(il),sp_efluct2(il), &
- sp_norm(il),sp_kinetic(il),sp_energy(il), &
- sp_orbital(:,il),sp_spin(:,il)
+ IF(il==npmin(2)) THEN
+ WRITE(*,'(A)') ' Proton Single Particle States:',header
+ END IF
+ WRITE(*,'(1X,I3,F4.0,F8.5,2F9.5,F9.6,F8.3,F10.3,3F8.3,3F7.3)') &
+ il,sp_parity(il),wocc(il),sp_efluct1(il),sp_efluct2(il), &
+ sp_norm(il),sp_kinetic(il),sp_energy(il), &
+ sp_orbital(:,il),sp_spin(:,il)
ENDDO
CALL moment_print
END SUBROUTINE sinfo
@@ -388,73 +388,73 @@ CONTAINS
!*************************************************************************
nst=0
DO iq=1,2
- IF(iq==1) THEN
- nps=npsi(1)
- ELSE
- nps=npsi(2)
- ENDIF
- ka_loop: DO ka=0,nps
- DO k=0,nps
- DO j=0,nps
- DO i=0,nps
- IF(ka==i+j+k) THEN
- DO is=1,2
- nst=nst+1
- IF(nst>nps) EXIT ka_loop
- nshell(1,nst)=i
- nshell(2,nst)=j
- nshell(3,nst)=k
+ IF(iq==1) THEN
+ nps=npsi(1)
+ ELSE
+ nps=npsi(2)
+ ENDIF
+ ka_loop: DO ka=0,nps
+ DO k=0,nps
+ DO j=0,nps
+ DO i=0,nps
+ IF(ka==i+j+k) THEN
+ DO is=1,2
+ nst=nst+1
+ IF(nst>nps) EXIT ka_loop
+ nshell(1,nst)=i
+ nshell(2,nst)=j
+ nshell(3,nst)=k
+ ENDDO
+ ENDIF
ENDDO
- ENDIF
- ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO ka_loop
- nst=nst-1
+ ENDDO ka_loop
+ nst=nst-1
END DO
DO iq=1,2
- nst=npmin(iq)
- DO iz=1,nz
- zz2=(z(iz)/radinz)**2
- DO ix=1,nx
- xx2=(x(ix)/radinx)**2
- DO iy=1,ny
- y2=(y(iy)/radiny)**2
- temp=xx2+y2+zz2
- psi(ix,iy,iz,1,nst)=EXP(-(temp))
+ nst=npmin(iq)
+ DO iz=1,nz
+ zz2=(z(iz)/radinz)**2
+ DO ix=1,nx
+ xx2=(x(ix)/radinx)**2
+ DO iy=1,ny
+ y2=(y(iy)/radiny)**2
+ temp=xx2+y2+zz2
+ psi(ix,iy,iz,1,nst)=EXP(-(temp))
+ ENDDO
ENDDO
- ENDDO
- ENDDO
- anorm=rpsnorm(psi(:,:,:,:,nst))
- psi(:,:,:,:,nst)=psi(:,:,:,:,nst)/SQRT(anorm)
- !*************************************************************************
- ! Higher states: lowest * polynomial
- !*************************************************************************
- DO nst=npmin(iq)+1,npsi(iq)
- is=MOD(nst-npmin(iq),2)+1
- DO iz=1,nz
- IF(nshell(3,nst)/=0) THEN
- zz=z(iz)**nshell(3,nst)
- ELSE
- zz=1.0D0
- ENDIF
- DO iy=1,ny
- IF(nshell(2,nst)/=0) THEN
- yy=y(iy)**nshell(2,nst)
- ELSE
- yy=1.0D0
- ENDIF
- DO ix=1,nx
- IF(nshell(1,nst)/=0) THEN
- xx=x(ix)**nshell(1,nst)
- ELSE
- xx=1.0D0
- ENDIF
- psi(ix,iy,iz,is,nst)=psi(ix,iy,iz,1,npmin(iq))*xx*yy*zz
- ENDDO
+ ENDDO
+ anorm=rpsnorm(psi(:,:,:,:,nst))
+ psi(:,:,:,:,nst)=psi(:,:,:,:,nst)/SQRT(anorm)
+ !*************************************************************************
+ ! Higher states: lowest * polynomial
+ !*************************************************************************
+ DO nst=npmin(iq)+1,npsi(iq)
+ is=MOD(nst-npmin(iq),2)+1
+ DO iz=1,nz
+ IF(nshell(3,nst)/=0) THEN
+ zz=z(iz)**nshell(3,nst)
+ ELSE
+ zz=1.0D0
+ ENDIF
+ DO iy=1,ny
+ IF(nshell(2,nst)/=0) THEN
+ yy=y(iy)**nshell(2,nst)
+ ELSE
+ yy=1.0D0
+ ENDIF
+ DO ix=1,nx
+ IF(nshell(1,nst)/=0) THEN
+ xx=x(ix)**nshell(1,nst)
+ ELSE
+ xx=1.0D0
+ ENDIF
+ psi(ix,iy,iz,is,nst)=psi(ix,iy,iz,1,npmin(iq))*xx*yy*zz
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- END DO
+ END DO
END DO
WRITE(*,*) "Harmonic oscillators widths (x-y-z):"
WRITE(*,"(3F12.4)") radinx,radiny,radinz
diff --git a/Code/twobody.f90 b/Code/twobody.f90
index 2e0576ecb58d1c682a3f46d2a4ae566c1838344b..1f17b906f35e99f9d5e3f4cfa0c027bcdef9372c 100644
--- a/Code/twobody.f90
+++ b/Code/twobody.f90
@@ -38,55 +38,55 @@ CONTAINS
tetold=ATAN(slold)
! *** Step 2
Iteration: DO itcm=1,10
- ! substep 1
- IF(itcm>1) THEN
- centx=centerx
- centz=centerz
- slope=(centz(2)-centz(1))/(centx(2)-centx(1))
- bb=centz(1)-slope*centx(1)
- ELSE ! first iteration
- centx=xold
- centz=zold
- CALL getslope
- bb=cmtot(3)-slope*cmtot(1)
- ENDIF
- ! substep 2
- istwobody=divpoint()
- ! substep 3
- angle=ATAN(slope)
- slopev=dtan(angle+pi/2.0D0)
- bb=zmin-slopev*xmin
- ! substep 4
- centerx=0.D0; centerz=0.D0
- mass=0.0D0
- charge=0.0D0
- DO iz=1,nz
- zz=z(iz)
- DO ix=1,nx
- xx=x(ix)
- diff=zz-slopev*xx-bb
- DO iy=1,ny
- vol=wxyz
- rhotot=rho(ix,iy,iz,1)+rho(ix,iy,iz,2)
- ifrag=1 ! 1 for left and 2 for right fragment
- IF(diff<0.0D0) THEN
- ifrag=1
- ELSE
- ifrag=2
- ENDIF
- mass(ifrag)=mass(ifrag)+vol*rhotot
- charge(ifrag)=charge(ifrag)+vol*rho(ix,iy,iz,2)
- centerx(ifrag)=centerx(ifrag)+vol*rhotot*xx
- centerz(ifrag)=centerz(ifrag)+vol*rhotot*zz
+ ! substep 1
+ IF(itcm>1) THEN
+ centx=centerx
+ centz=centerz
+ slope=(centz(2)-centz(1))/(centx(2)-centx(1))
+ bb=centz(1)-slope*centx(1)
+ ELSE ! first iteration
+ centx=xold
+ centz=zold
+ CALL getslope
+ bb=cmtot(3)-slope*cmtot(1)
+ ENDIF
+ ! substep 2
+ istwobody=divpoint()
+ ! substep 3
+ angle=ATAN(slope)
+ slopev=dtan(angle+pi/2.0D0)
+ bb=zmin-slopev*xmin
+ ! substep 4
+ centerx=0.D0; centerz=0.D0
+ mass=0.0D0
+ charge=0.0D0
+ DO iz=1,nz
+ zz=z(iz)
+ DO ix=1,nx
+ xx=x(ix)
+ diff=zz-slopev*xx-bb
+ DO iy=1,ny
+ vol=wxyz
+ rhotot=rho(ix,iy,iz,1)+rho(ix,iy,iz,2)
+ ifrag=1 ! 1 for left and 2 for right fragment
+ IF(diff<0.0D0) THEN
+ ifrag=1
+ ELSE
+ ifrag=2
+ ENDIF
+ mass(ifrag)=mass(ifrag)+vol*rhotot
+ charge(ifrag)=charge(ifrag)+vol*rho(ix,iy,iz,2)
+ centerx(ifrag)=centerx(ifrag)+vol*rhotot*xx
+ centerz(ifrag)=centerz(ifrag)+vol*rhotot*zz
+ ENDDO
ENDDO
- ENDDO
- ENDDO
- centerx=centerx/mass
- centerz=centerz/mass
- ! substep 5
- IF(.NOT.istwobody) EXIT
- IF(MAX(MAXVAL(ABS(centx-centerx)),MAXVAL(ABS(centz-centerz))) &
- <1.0d-05) EXIT
+ ENDDO
+ centerx=centerx/mass
+ centerz=centerz/mass
+ ! substep 5
+ IF(.NOT.istwobody) EXIT
+ IF(MAX(MAXVAL(ABS(centx-centerx)),MAXVAL(ABS(centz-centerz))) &
+ <1.0d-05) EXIT
ENDDO Iteration
!
! *** Step 3
@@ -111,19 +111,19 @@ CONTAINS
(ratio*(charge(1)*charge(2)*e2)**2))
ttt=xlf**2/(ratio*charge(1)*charge(2)*e2*roft)
IF(ABS(slope)>ABS(slold)) THEN
- tetc=pi-teti
- tetf=2.0D0*ACOS(1.0D0/epsf)-teti
+ tetc=pi-teti
+ tetf=2.0D0*ACOS(1.0D0/epsf)-teti
ELSE
- tetc=teti
- teti=pi-tetc
- tetf=ACOS(1.0D0/epsf)-ACOS(MAX(1.0D0/epsf *(1.0D0+ &
- ttt),1.D0))
+ tetc=teti
+ teti=pi-tetc
+ tetf=ACOS(1.0D0/epsf)-ACOS(MAX(1.0D0/epsf *(1.0D0+ &
+ ttt),1.D0))
ENDIF
tets=tetc-tetf
IF(ABS(slope) <ABS(slold)) THEN
- temp=mass(1); mass(1)=mass(2); mass(2)=temp
- temp=charge(1); charge(1)=charge(2); charge(2)=temp
- temp=tke2body(1); tke2body(1)=tke2body(2); tke2body(2)=temp
+ temp=mass(1); mass(1)=mass(2); mass(2)=temp
+ temp=charge(1); charge(1)=charge(2); charge(2)=temp
+ temp=tke2body(1); tke2body(1)=tke2body(2); tke2body(2)=temp
ENDIF
END SUBROUTINE twobody_case
!***************************************************************
@@ -132,20 +132,20 @@ CONTAINS
INTEGER :: ix,iy,iz
q2=0.D0
DO iz=1,nz
- zz=z(iz)-cmtot(3)
- DO iy=1,ny
- yy=y(iy)-cmtot(2)
- DO ix=1,nx
- xx=x(ix)-cmtot(1)
- vol=wxyz*(rho(ix,iy,iz,1)+rho(ix,iy,iz,2))
- q2(1,1)=q2(1,1)+(xx*xx+xx*xx-yy*yy-zz*zz)*vol
- q2(2,2)=q2(2,2)+(yy*yy+yy*yy-xx*xx-zz*zz)*vol
- q2(3,3)=q2(3,3)+(zz*zz+zz*zz-xx*xx-yy*yy)*vol
- q2(1,2)=q2(1,2)+3.D0*xx*yy*vol
- q2(1,3)=q2(1,3)+3.D0*xx*zz*vol
- q2(2,3)=q2(2,3)+3.D0*yy*zz*vol
- ENDDO
- ENDDO
+ zz=z(iz)-cmtot(3)
+ DO iy=1,ny
+ yy=y(iy)-cmtot(2)
+ DO ix=1,nx
+ xx=x(ix)-cmtot(1)
+ vol=wxyz*(rho(ix,iy,iz,1)+rho(ix,iy,iz,2))
+ q2(1,1)=q2(1,1)+(xx*xx+xx*xx-yy*yy-zz*zz)*vol
+ q2(2,2)=q2(2,2)+(yy*yy+yy*yy-xx*xx-zz*zz)*vol
+ q2(3,3)=q2(3,3)+(zz*zz+zz*zz-xx*xx-yy*yy)*vol
+ q2(1,2)=q2(1,2)+3.D0*xx*yy*vol
+ q2(1,3)=q2(1,3)+3.D0*xx*zz*vol
+ q2(2,3)=q2(2,3)+3.D0*yy*zz*vol
+ ENDDO
+ ENDDO
ENDDO
q2(2,1)=q2(1,2)
q2(3,1)=q2(1,3)
@@ -153,9 +153,9 @@ CONTAINS
denom=0.5D0*(q2(1,1)-q2(3,3)+SQRT((q2(1,1) &
-q2(3,3))**2+4.D0*q2(3,1)**2))
IF(ABS(denom)<1.D-4) THEN
- slope=100.D0
+ slope=100.D0
ELSE
- slope=q2(1,3)/denom
+ slope=q2(1,3)/denom
ENDIF
END SUBROUTINE getslope
!***************************************************************
@@ -173,54 +173,54 @@ CONTAINS
deltaz=0.5*dz
! Loop 1
IF(slope>=0.D0) THEN
- iz1=1; iz2=nz; idz=1
+ iz1=1; iz2=nz; idz=1
ELSE
- iz1=nz; iz2=1; idz=-1
+ iz1=nz; iz2=1; idz=-1
END IF
il=0
DO ix=1,nx
- xx=x(ix)
- DO iz=iz1,iz2,idz
- zz=z(iz)
- online=.FALSE.
- IF(ABS(slope)<=1.D0) THEN
- online=ABS(zz-slope*xx-bb)<=deltaz
- ELSE
- online=ABS(xx-(zz-bb)/slope)<=deltax
- END IF
- IF(online) THEN
- il=il+1
- IF(il>nx+ny+nz) THEN
- IF(wflag) WRITE(*,*) ' Increase dimensioning in function divpoint'
- STOP
+ xx=x(ix)
+ DO iz=iz1,iz2,idz
+ zz=z(iz)
+ online=.FALSE.
+ IF(ABS(slope)<=1.D0) THEN
+ online=ABS(zz-slope*xx-bb)<=deltaz
+ ELSE
+ online=ABS(xx-(zz-bb)/slope)<=deltax
+ END IF
+ IF(online) THEN
+ il=il+1
+ IF(il>nx+ny+nz) THEN
+ IF(wflag) WRITE(*,*) ' Increase dimensioning in function divpoint'
+ STOP
+ END IF
+ ixl(il)=ix; izl(il)=iz
+ rhol(il)=rho(ix,iyy,iz,1)+rho(ix,iyy,iz,2)
END IF
- ixl(il)=ix; izl(il)=iz
- rhol(il)=rho(ix,iyy,iz,1)+rho(ix,iyy,iz,2)
- END IF
- END DO
+ END DO
END DO
! Loop 2
nf=0; in_vacuum=.TRUE.
DO i=1,il
- IF(rhol(i)>vacuum) THEN
- IF(in_vacuum) THEN
- in_vacuum=.FALSE.
- nf=nf+1
- IF(nf==2) n2=i
- IF(nf>2) EXIT
- END IF
- ELSE
- IF(.NOT.in_vacuum) THEN
- IF(nf==1) n1=MAX(1,i-1)
- in_vacuum=.TRUE.
- END IF
- END IF
+ IF(rhol(i)>vacuum) THEN
+ IF(in_vacuum) THEN
+ in_vacuum=.FALSE.
+ nf=nf+1
+ IF(nf==2) n2=i
+ IF(nf>2) EXIT
+ END IF
+ ELSE
+ IF(.NOT.in_vacuum) THEN
+ IF(nf==1) n1=MAX(1,i-1)
+ in_vacuum=.TRUE.
+ END IF
+ END IF
END DO
two=nf==2
! final processing
IF(two) THEN
- xmin=0.5D0*(x(ixl(n1))+x(ixl(n2)))
- zmin=0.5D0*(z(izl(n1))+z(izl(n2)))
+ xmin=0.5D0*(x(ixl(n1))+x(ixl(n2)))
+ zmin=0.5D0*(z(izl(n1))+z(izl(n2)))
END IF
divpoint=two
END FUNCTION divpoint
@@ -240,20 +240,20 @@ CONTAINS
' tets(�):',tets*180.0D0/pi,' Vcoul:',xcoul, &
' Vcent:',xcent
IF(ABS(slope)>ABS(slold)) THEN
- centx=centerx; centz=centerz
+ centx=centerx; centz=centerz
ELSE
- centx(1)=centerx(2); centz(1)=centerz(2)
- centx(2)=centerx(1); centz(2)=centerz(1)
+ centx(1)=centerx(2); centz(1)=centerz(2)
+ centx(2)=centerx(1); centz(2)=centerz(1)
ENDIF
!
IF(wflag) THEN
- WRITE(*,'(/a/a)') ' Collision kinematics:', &
- ' Side Mass Charge <x> &
- &<y> tke'
- DO i=1,2
- WRITE(*,'(A,2F12.4,1P,3E12.4)') name(i),mass(i),charge(i), &
- centx(i),centz(i),tke2body(i)
- ENDDO
+ WRITE(*,'(/a/a)') ' Collision kinematics:', &
+ ' Side Mass Charge <x> &
+ &<y> tke'
+ DO i=1,2
+ WRITE(*,'(A,2F12.4,1P,3E12.4)') name(i),mass(i),charge(i), &
+ centx(i),centz(i),tke2body(i)
+ ENDDO
ENDIF
END SUBROUTINE twobody_print
END MODULE Twobody
diff --git a/Code/user_sample.f90 b/Code/user_sample.f90
index d6571d7945ca9adc2a4700a342872cf79706707a..8b6debda1064b96127e1146fb40f90b68d0422da 100644
--- a/Code/user_sample.f90
+++ b/Code/user_sample.f90
@@ -14,15 +14,15 @@ CONTAINS
,d,' and radii ',r
wocc=1.D0
DO ic=1,3
- FORALL(ix=1:nx,iy=1:ny,iz=1:nz)
- psi(ix,iy,iz,1,ic)= &
- EXP(-((x(ix)-center(ic))**2+y(iy)**2+z(iz)**2)/(2.D0*r**2))
- END FORALL
- psi(:,:,:,2,ic)=0.D0
- psi(:,:,:,1,ic+3)=0.D0
- psi(:,:,:,2,ic+3)=psi(:,:,:,1,ic)
- psi(:,:,:,:,ic+6)=psi(:,:,:,:,ic)
- psi(:,:,:,:,ic+9)=psi(:,:,:,:,ic+3)
+ FORALL(ix=1:nx,iy=1:ny,iz=1:nz)
+ psi(ix,iy,iz,1,ic)= &
+ EXP(-((x(ix)-center(ic))**2+y(iy)**2+z(iz)**2)/(2.D0*r**2))
+ END FORALL
+ psi(:,:,:,2,ic)=0.D0
+ psi(:,:,:,1,ic+3)=0.D0
+ psi(:,:,:,2,ic+3)=psi(:,:,:,1,ic)
+ psi(:,:,:,:,ic+6)=psi(:,:,:,:,ic)
+ psi(:,:,:,:,ic+9)=psi(:,:,:,:,ic+3)
END DO
END SUBROUTINE init_user
END MODULE User