The ABINIT Users Mailing List ( CLOSED )

[NISHIMATSU Takeshi <t-nissie@imr.tohoku.ac.jp>]

Dear Dr. Rushchanskii,

Here is my quick hack.
z-axis only, never stops before ntime.
Mr. Hashimoto, my co-worker, did more hacks based on it.
I am very grateful to authors for GPLed ABINIT. Thank you!

diff -uw src_tests_4.3.3-Alpha-original/Src_9drive/brdmin.f 
src_tests_4.3.3-Alpha-perovskite/Src_9drive/brdmin.f
--- src_tests_4.3.3-Alpha-original/Src_9drive/brdmin.f  2004-02-18 
00:34:39.000000000 +0900
+++ src_tests_4.3.3-Alpha-perovskite/Src_9drive/brdmin.f        2004-06-08 
15:52:24.000000000 +0900
@@ -162,7 +162,7 @@
  
 !Local variables-------------------------------
  integer,parameter :: level=5
- integer :: 
counter,cycle_main,hess_ok,iapp,iatom,iatom1,iatom2,idim,idir,idir1
+ integer :: 
counter,cycle_main,hess_ok,iapp,iatom,iatom1,iatom2,idim,idir,idir1,jdim
  integer :: idir2,iexit,ii,ios,itime,ixfh,jatom,jdir,jj,ndim,ndim0,openexit
  integer :: optcell,option,prtvel,prtvol,routine
  real(dp) :: DDOT,diag,etotal_prev,favg,ucvol,ucvol0
@@ -171,8 +171,9 @@
  character*500 :: message
  real(dp) :: acell0(3),angle(3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
  real(dp) :: rprimd0(3,3)
- real(dp),allocatable :: fred_corrected(:,:),hessin(:,:),vin(:),vin_prev(:)
- real(dp),allocatable :: vout(:),vout_prev(:),xcart(:,:)
+ real(dp),allocatable :: 
fred_corrected(:,:),hessin(:,:),vin(:),vin_prev(:),n_vin(:),n_hessin(:)
+ real(dp),allocatable :: vout(:),vout_prev(:),xcart(:,:),modified_vout(:)
+ real(dp) :: u_z, amp, ratio, norm, prod, lambda, denominator
  
 !***************************************************************************
 !Beginning of executable session
@@ -195,7 +196,7 @@
  if(optcell==7 .or. optcell==8 .or. optcell==9)ndim=ndim+3
 
  allocate(fred_corrected(3,natom),xcart(3,natom))
- allocate(vin(ndim),vout(ndim),vin_prev(ndim),vout_prev(ndim))
+ 
allocate(vin(ndim),vout(ndim),vin_prev(ndim),vout_prev(ndim),modified_vout(ndim),n_vin(ndim),n_hessin(ndim))
  allocate(hessin(ndim,ndim))
 
 !Compute dimensional primitive translations rprimd, then metric tensor gmet
@@ -320,11 +321,75 @@
     etotal_prev=results_gs%etotal
     vin_prev(:)=vin(:)
  
+    !t-nissie start
+    write(message,'(a,f21.16)') 'previous acell(3)=', acell(3)
+    call wrtout(ab_out,message,'COLL')
+    call perovskite_remove_translation(vin)
+    call perovskite_displacement_norm(vin,norm)
+    if(itime==1)then
+     u_z = acell(3) * norm
+     write(message,'(a,f21.16)') 'u_z=', u_z
+     call wrtout(ab_out,message,'COLL')
+    end if
+    write(message,'(a,7f14.10)') 'BEFORE vin, acell=', vin(3), vin(6), 
vin(9), vin(12), vin(15), acell(1), acell(3)
+    call wrtout(ab_out,message,'COLL')
+    write(message,'(a,7f14.10)') '             vout=', vout(3), vout(6), 
vout(9), vout(12), vout(15), vout(16), vout(18)
+    call wrtout(ab_out,message,'COLL')
+
+    n_vin=0.0d0
+    n_vin( 3) =  vin( 3)       /norm
+    n_vin( 6) = (vin( 6)-0.5d0)/norm
+    n_vin( 9) = (vin( 9)-0.5d0)/norm
+    n_vin(12) = (vin(12)-0.5d0)/norm
+    n_vin(15) =  vin(15)       /norm
+    write(message,'(a,5f14.10)') '            n_vin=', n_vin(3), n_vin(6), 
n_vin(9), n_vin(12), n_vin(15)
+    call wrtout(ab_out,message,'COLL')
+    n_hessin=0.0d0
+    do jdim=1,ndim
+       do idim=1,ndim
+          n_hessin(jdim)=n_hessin(jdim)+n_vin(idim)*hessin(idim,jdim)
+       enddo
+    enddo
+    lambda=DOT_PRODUCT(n_hessin(:),vout(:))
+    denominator=DOT_PRODUCT(n_hessin(:),n_vin(:))
+    lambda=lambda/denominator
+    write(message,'(a,f21.16)') 'lambda=', lambda
+    call wrtout(ab_out,message,'COLL')
+    modified_vout(:) = vout(:)
+    modified_vout( 3) = modified_vout( 3) - lambda*n_vin( 3)
+    modified_vout( 6) = modified_vout( 6) - lambda*n_vin( 6)
+    modified_vout( 9) = modified_vout( 9) - lambda*n_vin( 9)
+    modified_vout(12) = modified_vout(12) - lambda*n_vin(12)
+    modified_vout(15) = modified_vout(15) - lambda*n_vin(15)
+    write(message,'(a,7f14.10)') '    modified_vout=', modified_vout(3), 
modified_vout(6), modified_vout(9), modified_vout(12), modified_vout(15), 
modified_vout(16), modified_vout(18)
+    call wrtout(ab_out,message,'COLL')
+    !t-nissie end
+
 !   New atomic cartesian coordinates are obtained from vin, hessin and vout
     do idim=1,ndim
-     vin(:)=vin(:)-hessin(:,idim)*vout(idim)
+     vin(:)=vin(:)-hessin(:,idim)*modified_vout(idim)
     enddo
  
+    !t-nissie start   THIS SHOULD BE IMPLEMENTED IN Src_5common/xfpack.f
+    write(message,'(a,7f14.10)') 'MID    vin, acell=', vin(3), vin(6), 
vin(9), vin(12), vin(15), acell0(1)*vin(16), acell0(3)*vin(18)
+    call wrtout(ab_out,message,'COLL')
+    call perovskite_remove_translation(vin)
+    call perovskite_displacement_norm(vin,norm)
+    amp = acell0(3)*vin(18)*norm
+    write(message,'(a,f21.16)') 'amp=', amp
+    call wrtout(ab_out,message,'COLL')
+    ratio = u_z/amp
+    write(message,'(a,f21.16)') 'ratio=', ratio
+    call wrtout(ab_out,message,'COLL')
+    vin( 3) =         (vin( 3)      )*ratio
+    vin( 6) = 0.5d0 + (vin( 6)-0.5d0)*ratio
+    vin( 9) = 0.5d0 + (vin( 9)-0.5d0)*ratio
+    vin(12) = 0.5d0 + (vin(12)-0.5d0)*ratio
+    vin(15) =         (vin(15)      )*ratio
+    write(message,'(a,7f14.10)') 'AFTER  vin, acell=', vin(3), vin(6), 
vin(9), vin(12), vin(15), acell0(1)*vin(16), acell0(3)*vin(18)
+    call wrtout(ab_out,message,'COLL')
+    !t-nissie end   THIS SHOULD BE IMPLEMENTED IN Src_5common/xfpack.f
+
 !   Previous atomic forces
     vout_prev(:)=vout(:)
 
@@ -595,3 +660,27 @@
 
  end subroutine
 !!***
+
+subroutine perovskite_remove_translation(vin)
+  use defs_basis
+  implicit none
+  character*500 :: message
+  real(dp) :: vin(*), translation
+  translation = (vin(3) + (vin(6)-0.5d0) + (vin(9)-0.5d0) + (vin(12)-0.5d0) 
+ vin(15))/5.0d0
+  write(message,'(a,f21.16,a)') 'translation ', translation, ' is removed.'
+  call wrtout(ab_out,message,'COLL')
+  vin( 3) = vin( 3)-translation
+  vin( 6) = vin( 6)-translation
+  vin( 9) = vin( 9)-translation
+  vin(12) = vin(12)-translation
+  vin(15) = vin(15)-translation
+end subroutine perovskite_remove_translation
+
+subroutine perovskite_displacement_norm(vin,norm)
+  use defs_basis
+  implicit none
+  real(dp) :: vin(*),norm
+  norm = sqrt(vin(3)**2 + (vin( 6)-0.5d0)**2 &
+       &                + (vin( 9)-0.5d0)**2 &
+       &                + (vin(12)-0.5d0)**2 + vin(15)**2)
+end subroutine perovskite_displacement_norm


# Input (PbTiO3):
optcell 2
ionmov 2
ntime 30

prtwf 0

ecut 33.0
ecutsm 3.0
dilatmx 1.1
ixc 3

acell 3.9000 3.9000 3.9000 angstrom

rprim  1 0 0
       0 1 0
       0 0 1

#Definition of the atom types
ntypat 3
znucl 82 22 8

#Definition of the atoms
natom 5
typat 1*1 1*2 3*3
xred
  0.0 0.0 0.041324
  0.5 0.5 0.524248
  0.0 0.5 0.473788
  0.5 0.0 0.473788
  0.5 0.5 -0.013148

#Definition of the k-point grid
kptopt 1
ngkpt  8 8 8
nshiftk 1
shiftk 0.5 0.5 0.5

#Definition of the SCF procedure
nstep 35
toldff 5.0d-7     # Use force!
diemac 12.0



Sincerely,
-- 
 love && peace && free_software
 NISHIMATSU Takeshi