root / src / CheckPeriodicBound.f90 @ 3
Historique | Voir | Annoter | Télécharger (2,54 ko)
| 1 |
SUBROUTINE CheckPeriodicBound |
|---|---|
| 2 |
! This subroutine is only for periodic systems |
| 3 |
! it checks that the atoms do not change their |
| 4 |
! positions from one side to the other during the simulations. |
| 5 |
! This would ruin the interpolation. |
| 6 |
|
| 7 |
|
| 8 |
Use Path_module |
| 9 |
Use Io_module |
| 10 |
|
| 11 |
IMPLICIT NONE |
| 12 |
|
| 13 |
|
| 14 |
INTEGER(KINT) :: IGeom, I,J,K |
| 15 |
REAL(KREAL), ALLOCATABLE :: GeomRef(:,:), GeomTmp(:,:) ! (3,Nat) |
| 16 |
CHARACTER(SCHARS) :: Line,Line2 |
| 17 |
|
| 18 |
REAL(KREAL) :: Latrloc(3,3),Bloc(3),Xt,Yt,Zt |
| 19 |
LOGICAL :: Debug |
| 20 |
|
| 21 |
INTERFACE |
| 22 |
function valid(string) result (isValid) |
| 23 |
CHARACTER(*), intent(in) :: string |
| 24 |
logical :: isValid |
| 25 |
END function VALID |
| 26 |
END INTERFACE |
| 27 |
|
| 28 |
debug=valid('checkperiodicbound')
|
| 29 |
if (debug) WRITE(*,*) "========================= Entering CheckPeriodicBound ===========" |
| 30 |
|
| 31 |
ALLOCATE(GeomRef(3,Nat), GeomTmp(3,Nat)) |
| 32 |
|
| 33 |
V_direct=adjustl(Vasp_direct) |
| 34 |
Call upcase(V_direct) |
| 35 |
|
| 36 |
! The test is easy in DIRECT space, so we will convert the geometries |
| 37 |
! into this space. |
| 38 |
! We might optimize this as Read_geom converts geometries into NORMAL |
| 39 |
! coordinates... so it might be that we are converting them back to direct here... |
| 40 |
|
| 41 |
Latrloc(1:3,1)=Lat_a |
| 42 |
Latrloc(1:3,2)=Lat_b |
| 43 |
Latrloc(1:3,3)=Lat_c |
| 44 |
Bloc=1. |
| 45 |
CALL Gaussj(Latrloc,3,3,Bloc,1,1) |
| 46 |
|
| 47 |
! We put first geometry as the reference |
| 48 |
GeomRef=0. |
| 49 |
DO I=1,Nat |
| 50 |
DO K=1,3 |
| 51 |
DO J=1,3 |
| 52 |
GeomRef(K,I)=GeomRef(K,I)+XyzGeomI(1,J,I)*LatrLoc(K,J) |
| 53 |
END DO |
| 54 |
END DO |
| 55 |
END DO |
| 56 |
|
| 57 |
DO IGeom=2,NGeomI |
| 58 |
GeomTmp=0. |
| 59 |
DO I=1,Nat |
| 60 |
DO K=1,3 |
| 61 |
DO J=1,3 |
| 62 |
GeomTmp(K,I)=GeomTmp(K,I)+XyzGeomI(IGeom,J,I)*LatrLoc(K,J) |
| 63 |
END DO |
| 64 |
if (abs(GeomTmp(K,I)-GeomRef(K,I)).GE.BoxTol) THEN |
| 65 |
If (debug) & |
| 66 |
WRITE(*,'("Coord ",I1," has been changed for atom ",I5," by ",F4.0," ref:",F9.6," old:",F9.6," new:",F9.6)') &
|
| 67 |
K,I,Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)),& |
| 68 |
GeomRef(K,I),GeomTmp(K,I), & |
| 69 |
GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)) |
| 70 |
GeomTmp(K,I)=GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)) |
| 71 |
END IF |
| 72 |
END DO |
| 73 |
Xt=GeomTmp(1,I) |
| 74 |
Yt=GeomTmp(2,I) |
| 75 |
Zt=GeomTmp(3,I) |
| 76 |
XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1) |
| 77 |
XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2) |
| 78 |
XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3) |
| 79 |
END DO |
| 80 |
END DO |
| 81 |
|
| 82 |
if (debug) WRITE(*,*) "========================= Exiting CheckPeriodicBound ===========" |
| 83 |
END SUBROUTINE CheckPeriodicBound |
| 84 |
|
| 85 |
|
| 86 |
|
| 87 |
|
| 88 |
|
| 89 |
|