diff --git a/Code/Changes.txt b/Code/Changes.txt
index cf57e2ab9b2437655199fb32ac39882a30c5daee..4262b8c2f579ddc1dee852ee87776f3ee4dd58ef 100644
--- a/Code/Changes.txt
+++ b/Code/Changes.txt
@@ -1 +1,37 @@
-Corrected Sly4d in forces.data.
+Changes relative to published CPC version from April 2014:
+----------------------------------------------------------
+
+
+Step 1.0:
+
+- Fortran90 standard for indentation.
+
+- dynamic.f90: Enforce two-body analysis and writeing of restart file
+ at end of run, eben if not multiple of 'mrest'.
+
+- forces.data: SLy4d corrected
+
+- main3d.f90: New variable "nofsave" to maintain two-body structure
+ after restart.
+
+
+
+Step 2.0:
+
+- energies.f90: New variables to store previous energy fluctuationas.
+
+- main3d.f90: No new ortho-normalization after restart.
+
+- meanfield.f90: 'workden', 'workvec', 'pwwk', 'pswk2' made ALLOCATABLE.
+ Stabilized treatment of l*s potential. Better but
+ slower. Old version kept in 'meanfield-old.f90'.
+
+- static.f90: Accelerated iteration by flexible step size.
+ 'npsi' always set to even for consistent pairing.
+
+
+
+
+
+
+
diff --git a/Code/energies.f90 b/Code/energies.f90
index a045175f6342d0103b6fdae1212e8c50b93e14b3..404fadd02a5f4d99643e4cdedfb265fd2fa187eb 100644
--- a/Code/energies.f90
+++ b/Code/energies.f90
@@ -18,9 +18,12 @@ MODULE Energies
REAL(db) :: ecorc ! Slater & Koopman exchange
REAL(db) :: ehfint ! integrated total energy
REAL(db) :: efluct1 ! energyfluctuation h**2
+ REAL(db) :: efluct1prev
REAL(db) :: efluct2 ! fluctuation h*efluct
+ REAL(db) :: efluct2prev
REAL(db) :: tke ! kinetic energy summed
REAL(db) :: ehf ! Hartree-Fock energy from s.p. levels
+ REAL(db) :: ehfprev
REAL(db) :: e3corr ! rearrangement energy
REAL(db) :: orbital(3),spin(3),total_angmom(3)
CONTAINS
diff --git a/Code/main3d.f90 b/Code/main3d.f90
index bf825ac0acf9c003b5fda163564b1b98f914dc3e..c8c73b21143e37f487cbb416de8d158a806cdf97 100644
--- a/Code/main3d.f90
+++ b/Code/main3d.f90
@@ -114,7 +114,7 @@ PROGRAM tdhf3d
!********************************************************************
IF(nof>0) THEN
CALL read_fragments
- IF(.NOT.tmpi) THEN
+ IF(.NOT.(tmpi.OR.trestart)) THEN
CALL schmid
WRITE(*,*) 'Reorthogonalization complete'
END IF
diff --git a/Code/meanfield-old.f90 b/Code/meanfield-old.f90
new file mode 100644
index 0000000000000000000000000000000000000000..97eb15f3a24143004a6d84920905af3842704229
--- /dev/null
+++ b/Code/meanfield-old.f90
@@ -0,0 +1,208 @@
+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) :: workden(nx,ny,nz,2),workvec(nx,ny,nz,3,2)
+ INTEGER :: ix,iy,iz,ic,iq,icomp
+ ! 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
+ 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),DIMENSION(nx,ny,nz,2) :: pswk,pswk2
+ ! Step 1: non-derivative parts not involving spin
+ DO is=1,2
+ pout(:,:,:,is)=CMPLX(upot(:,:,:,iq)-eshift, &
+ -0.5D0*divaq(:,:,:,iq),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
+ pout(:,:,:,is)=pout(:,:,:,is) &
+ -CMPLX(dbmass(:,:,:,1,iq),aq(:,:,:,1,iq) &
+ -sigis*wlspot(:,:,:,2,iq),db)*pswk(:,:,:,is) &
+ -sigis*wlspot(:,:,:,3,iq)*pswk(:,:,:,ic) &
+ -bmass(:,:,:,iq)*pswk2(:,:,:,is)
+ ENDDO
+ ! 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
+ pout(:,:,:,is)=pout(:,:,:,is) &
+ -CMPLX(dbmass(:,:,:,2,iq),aq(:,:,:,2,iq) &
+ +sigis*wlspot(:,:,:,1,iq),db)*pswk(:,:,:,is) &
+ +CMPLX(0.D0,wlspot(:,:,:,3,iq),db)*pswk(:,:,:,ic) &
+ -bmass(:,:,:,iq)*pswk2(:,:,:,is)
+ ENDDO
+ ! 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
+ pout(:,:,:,is)=pout(:,:,:,is) &
+ -CMPLX(dbmass(:,:,:,3,iq),aq(:,:,:,3,iq),db)*pswk(:,:,:,is) &
+ +CMPLX(sigis*wlspot(:,:,:,1,iq),-wlspot(:,:,:,2,iq),db)* &
+ pswk(:,:,:,ic)-bmass(:,:,:,iq)*pswk2(:,:,:,is)
+ ENDDO
+ END SUBROUTINE hpsi
+ !***********************************************************************
+END Module Meanfield
diff --git a/Code/meanfield.f90 b/Code/meanfield.f90
index 97eb15f3a24143004a6d84920905af3842704229..462b626955f53365d8aa1a21214aa288df96d12d 100644
--- a/Code/meanfield.f90
+++ b/Code/meanfield.f90
@@ -26,8 +26,9 @@ CONTAINS
USE Trivial, ONLY: rmulx,rmuly,rmulz
REAL(db),PARAMETER :: epsilon=1.0d-25
REAL(db) :: rotspp,rotspn
- REAL(db) :: workden(nx,ny,nz,2),workvec(nx,ny,nz,3,2)
+ 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
@@ -130,6 +131,7 @@ CONTAINS
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)
@@ -143,11 +145,12 @@ CONTAINS
INTENT(OUT) :: pout
INTEGER :: is,ic
REAL(db) :: sigis
- COMPLEX(db),DIMENSION(nx,ny,nz,2) :: pswk,pswk2
+ 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, &
- -0.5D0*divaq(:,:,:,iq),db)*pinn(:,:,:,is)
+ -0D0,db)*pinn(:,:,:,is)
ENDDO
! Step 2: the spin-current coupling
pout(:,:,:,1)=pout(:,:,:,1) &
@@ -165,13 +168,25 @@ CONTAINS
ENDIF
DO is=1,2
ic=3-is
- sigis=3-2*is
+ sigis=(3-2*is)*0.5D0
pout(:,:,:,is)=pout(:,:,:,is) &
- -CMPLX(dbmass(:,:,:,1,iq),aq(:,:,:,1,iq) &
+ -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)
@@ -181,13 +196,25 @@ CONTAINS
ENDIF
DO is=1,2
ic=3-is
- sigis=3-2*is
+ sigis=(3-2*is)*0.5D0
pout(:,:,:,is)=pout(:,:,:,is) &
- -CMPLX(dbmass(:,:,:,2,iq),aq(:,:,:,2,iq) &
+ -CMPLX(dbmass(:,:,:,2,iq),0.5D0*aq(:,:,:,2,iq) &
+sigis*wlspot(:,:,:,1,iq),db)*pswk(:,:,:,is) &
- +CMPLX(0.D0,wlspot(:,:,:,3,iq),db)*pswk(:,:,:,ic) &
+ +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)
@@ -197,12 +224,24 @@ CONTAINS
ENDIF
DO is=1,2
ic=3-is
- sigis=3-2*is
+ sigis=(3-2*is)*0.5D0
pout(:,:,:,is)=pout(:,:,:,is) &
- -CMPLX(dbmass(:,:,:,3,iq),aq(:,:,:,3,iq),db)*pswk(:,:,:,is) &
- +CMPLX(sigis*wlspot(:,:,:,1,iq),-wlspot(:,:,:,2,iq),db)* &
+ -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/static.f90 b/Code/static.f90
index f5f7230f414f2c7c5fe14dcef610451948e95ee5..6b6df56d85c1370f71b88022bd65eb8b08e38899 100644
--- a/Code/static.f90
+++ b/Code/static.f90
@@ -12,14 +12,15 @@ MODULE Static
IMPLICIT NONE
LOGICAL :: tdiag=.FALSE.
LOGICAL :: tlarge=.FALSE.
+ LOGICAL :: tvaryx_0=.FALSE.
INTEGER :: maxiter
REAL(db) :: radinx,radiny,radinz, &
- serr,delesum,x0dmp=0.2D0,e0dmp=100.D0
+ serr,delesum,x0dmp=0.2D0,e0dmp=100.D0,x0dmpmin=0.2d0
CONTAINS
!*************************************************************************
SUBROUTINE getin_static
NAMELIST/static/ tdiag,tlarge,maxiter, &
- radinx,radiny,radinz,serr,x0dmp,e0dmp,nneut,nprot,npsi
+ radinx,radiny,radinz,serr,x0dmp,e0dmp,nneut,nprot,npsi,tvaryx_0
npsi=0
READ(5,static)
IF(nof<=0) THEN
@@ -28,6 +29,7 @@ CONTAINS
npsi(1)=nneut
ELSE
npsi(1)=NINT(nneut+1.65*FLOAT(nneut)**0.666667D0)
+ IF(MOD(npsi(1),2)/=0) npsi(1)=npsi(1)+1
ENDIF
ENDIF
IF(npsi(2)==0) THEN
@@ -35,6 +37,7 @@ CONTAINS
npsi(2)=nprot
ELSE
npsi(2)=NINT(nprot+1.65*FLOAT(nprot)**0.666667D0)
+ IF(MOD(npsi(2),2)/=0) npsi(2)=npsi(2)+1
ENDIF
ENDIF
IF(nneut>npsi(1).OR.nprot>npsi(2)) &
@@ -42,6 +45,8 @@ CONTAINS
nstmax=npsi(1)+npsi(2)
charge_number=nprot
mass_number=nneut+nprot
+ x0dmpmin=x0dmp
+ WRITE(*,*) "x0dmpmin=", x0dmpmin
END IF
END SUBROUTINE getin_static
!*************************************************************************
@@ -122,9 +127,11 @@ CONTAINS
! produce and print detailed information
CALL sp_properties
CALL sinfo
+ !set x0dmp to 3* its value to get faster convergence
+ IF(tvaryx_0) x0dmp=3.0d0*x0dmp
! step 4: start static iteration loop
Iteration: DO iter=firstiter,maxiter
- WRITE(*,'(a,i6)') ' Static Iteration No.',iter
+ WRITE(*,'(a,i6,a,F12.4)') ' Static Iteration No.',iter,' x0dmp= ',x0dmp
! Step 5: gradient step
delesum=0.0D0
sumflu=0.0D0
@@ -171,9 +178,12 @@ CONTAINS
ENDIF
CALL skyrme
! calculate and print information
- IF(mprint>0.AND.MOD(iter,mprint)==0) THEN
+ IF(mprint>0.AND.MOD(iter,mprint)==0) THEN
CALL sp_properties
CALL sinfo
+ ELSEIF(tvaryx_0) THEN
+ CALL sp_properties
+ CALL sum_energy
ENDIF
! Step 9: check for convergence, saving wave functions
IF(sumflu/nstmax<serr.AND.iter>1) THEN
@@ -183,6 +193,19 @@ CONTAINS
IF(MOD(iter,mrest)==0) THEN
CALL write_wavefunctions
ENDIF
+ ! Step 10: update step size for the next iteration
+ IF(tvaryx_0) THEN
+ IF(ehf<ehfprev .OR. efluct1<(efluct1prev*(1.0d0-1.0d-5)) &
+ .OR. efluct2<(efluct2prev*(1.0d0-1.0d-5))) THEN
+ x0dmp=x0dmp*1.005
+ ELSE
+ x0dmp=x0dmp*0.8
+ END IF
+ IF(x0dmp<x0dmpmin) x0dmp=x0dmpmin
+ efluct1prev=efluct1
+ efluct2prev=efluct2
+ ehfprev=ehf
+ END IF
END DO Iteration
IF(tdiag) DEALLOCATE(hmatr)
END SUBROUTINE statichf
@@ -216,7 +239,7 @@ CONTAINS
ENDDO
IF(tdiag) hmatr(nst,nst)=hmatr(nst,nst)+spe
! Step 3: calculate fluctuation, i.e. <h*h> and |h|**2
- IF(mprint>0.AND.MOD(iter,mprint)==0) THEN
+ IF((mprint>0.AND.MOD(iter,mprint)==0).OR.tvaryx_0) THEN
CALL hpsi(iq,esf,ps1,ps2)
exph2=overlap(psin,ps2)
varh2=rpsnorm(ps1)
@@ -340,12 +363,12 @@ CONTAINS
OPEN(unit=scratch,file=spinfile, POSITION='APPEND')
WRITE(scratch,'(1x,i5,9F10.4)') iter,orbital,spin,total_angmom
CLOSE(unit=scratch)
- WRITE(*,'(/,A,I7,A/2(A,F12.4),A/(3(A,E12.5),A))') &
+ WRITE(*,'(/,A,I7,A/2(A,F12.4),A/(3(A,E12.5),A),3(A,E12.5))') &
' ***** Iteration ',iter,' *****',' Total energy: ',ehf,' MeV Total kinetic energy: ', &
tke,' MeV',' de/e: ',delesum,' h**2 fluct.: ',efluct1, &
' MeV, h*hfluc.: ',efluct2,' MeV', &
' MeV. Rearrangement E: ',e3corr,' MeV. Coul.Rearr.: ', &
- ecorc,' MeV'
+ ecorc,' MeV x0dmp: ',x0dmp
! detail printout
WRITE(*,'(/A)') ' Energies integrated from density functional:'
WRITE(*,'(4(A,1PE14.6),A/26X,3(A,1PE14.6),A)') &