## Revision 10

3 3
  ! This routines calculates the energy and the gradient of

4 4
  ! a molecule, using mopac

5 5

6
  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, ProgExe

6
  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, ProgExe, &

7
       FPBC, Lat_a, Lat_b, Lat_c, IPer

7 8
  use Io_module

8 9

9 10
  IMPLICIT NONE

......
92 93
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:)

93 94
     END IF

94 95
  END DO

96
! If we have periodic boundaries, we have to print them

97
  IF (FPBC) THEN

98
     IF (IPER>=1) THEN

99
        WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_a(1:3)

100
     END IF

101
     IF (IPER>=2) THEN

102
        WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_b(1:3)

103
     END IF

104
     IF (IPER>=3) THEN

105
        WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_c(1:3)

106
     END IF

107
  END IF

95 108

96 109
     WRITE(IOTMP,'(1X,A)') Trim(Mopac_EndGeom)

97 110

src/AlignPartial.f90 (revision 10)
13 13
  ! x2,y2,z2  : coordinates of the second geometry

14 14
  ! ListAt    : Logical, T if this atom should be considered in the RMSD calc and

15 15
  !             alignment.

16
  ! Debug     : True if additionnal printing is requiered

17 16
  !

18 17
  ! output parameters:

19 18
  ! U         : Real, (3x3) matrix for rotation

20
  ! Rmsd      : Read, the rmsd

21 19
  ! x2,y2,z2  : coordinates of the second geometry ALIGNED to the first one

22 20
  !-----------------------------------------------------------------------

23 21

13 13
       CHARACTER(*), intent(in) :: string

14 14
       logical                  :: isValid

15 15
     END function VALID

16

17
    SUBROUTINE die(routine, msg, file, line, unit)

18

19
      Use VarTypes

20
      Use io_module

21

22
      implicit none

23
      character(len=*), intent(in)           :: routine, msg

24
      character(len=*), intent(in), optional :: file

25
      integer(KINT), intent(in), optional      :: line, unit

26

27
    END SUBROUTINE die

28

29

16 30
  END INTERFACE

17 31

18 32

19
  CHARACTER(132) ::  Line,Line2

33
  CHARACTER(LCHARS) ::  Line,LineUp

20 34
  INTEGER(KINT) :: LineL, Idx, NTmp

35
  INTEGER(KINT) :: NatMopac

36
  REAL(KREAL) :: Lat(3,3)

21 37


22 38
  LOGICAL :: Debug

23 39

......
50 66
     LineL=len(Trim(Line))

51 67
     IF (Line(1:1)/="*") THEN

52 68
        IF (NTmp==0) THEN

53
           Line2=Line

54
           Call UpCase(Line2)

55
           Idx=Index(Line2,'GRADIENTS')

69
           LineUp=Line

70
           Call UpCase(LineUp)

71
           Idx=Index(LineUp,'GRADIENTS')

56 72
           If (Idx==0) Line=TRIM(Line) // " GRADIENTS"

57
           Idx=Index(Line2,'1SCF')

73
           Idx=Index(LineUp,'1SCF')

58 74
           If (Idx==0) Line=TRIM(Line) // " 1SCF"

59 75
        END IF

60 76
        current%Line=TRIM(Line)

......
77 93
!     END DO

78 94

79 95
! Now the geometry... that we just skip :)

96
! PFL 2013 Apr

97
! We take care that there is no Translation vectors...

98
! We also check that the number of atoms is ok

80 99
  IF (DEBUG) WRITE(*,*) "Reading Mopac Geometry"

81 100
  Mopac_EndGeom=""

82 101
  LineL=1

102
  NatMopac=0

103
  Lat=0.d0

104
  IPer=0

105
  FPBC=.FALSE.

83 106
  DO WHILE (LineL.NE.0)

84 107
     READ(IOIN,'(A)',END=989) Line

85 108
     Line=AdjustL(Line)

86 109
     LineL=len(Trim(Line))

87 110
     ! The last line might be either blank or filled with 0

88
     IF (Line(1:1)=="0") THEN

89
        LineL=0

90
        Mopac_EndGeom=Trim(Line)

111
     If (LineL>0) THEN

112
        SELECT CASE (Line(1:1))

113
          CASE ("0")

114
             LineL=0

115
             Mopac_EndGeom=Trim(Line)

116
          CASE("*")

117
             CurCom%Line=TRIM(LINE)

118
             ALLOCATE(CurCom%Next)

119
             CurCom => CurCom%Next

120
             NULLIFY(CurCom%Next)

121
          CASE DEFAULT

122
             LineUp=Line

123
             Call UpCase(LineUp)

124
             If (LineUp(1:2)=="TV") THEN

125
                FPBC=.TRUE.

126
                IPer=IPer+1

127
                If (Iper>3) THEN

128
                   Call Die("ReadInput Mopac","Iper>3",Unit=IOOUT)

129
                END IF

130
                NTmp=Index(LineUp," ")

131
                LineUp=LineUp(NTmp:)

132
                Read(LineUp,*) Lat(IPer,1:3)

133
             ELSE

134
                NatMopac=NatMopac+1

135
             END IF

136
         END SELECT

137
      END IF

138


139
   END DO

140

141
!  WRITE(*,*) "NatMopac,Nat:",NAtMopac,Nat

142
  IF (NatMopac/=Nat) Call Die("ReadInput_mopac","Nat read does not mat nat",Unit=IOOUT)

143
  IF (FPBC) THEN

144
     Lat_a(1:3)=Lat(1,1:3)

145
     Lat_b(1:3)=Lat(2,1:3)

146
     Lat_c(1:3)=Lat(3,1:3)

147
     If (IPer>=1) THEN

148
        kaBeg=-1

149
        kaEnd=1

91 150
     END IF

92
     IF (Line(1:1)=="*") THEN

93
        CurCom%Line=TRIM(LINE)

94
        ALLOCATE(CurCom%Next)

95
        CurCom => CurCom%Next

96
        NULLIFY(CurCom%Next)

151
     If (IPer>=2) THEN

152
        kbBeg=-1

153
        kbEnd=1

97 154
     END IF

98
  END DO

155
     If (IPer==3) THEN

156
        kcBeg=-1

157
        kcEnd=1

158
     END IF

159
     If (IPer>3) THEN

160
        Call Die("Readinput_mopac","Found too many Tv lines !",Unit=IOOUT)

161
     END IF

162
  END IF

99 163

100 164
! If we are here, there might be something else to read: Mopac_end

101 165

src/Io_module.f90 (revision 10)
10 10
  INTEGER(KINT) :: IOIN=11, IOOUT=12, IOCART=14

11 11
  INTEGER(KINT) :: IOGEOM=15, IODAT=16,IoGplot=17

12 12
  INTEGER(KINT), PARAMETER :: IOTMP=21,IOTMP2=22, IOTMP3=23

13
  INTEGER(KINT), PARAMETER :: IOERR=19

13
  INTEGER(KINT), PARAMETER :: IOERR=0

14 14
  CHARACTER(SCHARS) :: RunMode

15 15

16 16
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

13 13
       CHARACTER(*), intent(in) :: string

14 14
       logical                  :: isValid

15 15
     END function VALID

16

17

18
    SUBROUTINE die(routine, msg, file, line, unit)

19

20
      Use VarTypes

21
      Use io_module

22

23
      implicit none

24
      character(len=*), intent(in)           :: routine, msg

25
      character(len=*), intent(in), optional :: file

26
      integer(KINT), intent(in), optional      :: line, unit

27

28
    END SUBROUTINE die

29

30

16 31
  END INTERFACE

17 32

18 33

19
  CHARACTER(132) ::  Line

34
  CHARACTER(132) ::  Line,LineUp

20 35
  INTEGER(KINT) :: LineL, Idx, Iat

21
  INTEGER(KINT) :: I

36
  INTEGER(KINT) :: I,NTmp

37
  REAL(KREAL) :: Lat(3,3)

22 38


23 39
  LOGICAL :: Debug

24 40

......
86 102
     ALLOCATE(Gauss_Paste(NAt))

87 103
     LineL=1

88 104
     Iat=0

105
     IPer=0

106
     FPBC=.FALSE.

89 107
     Gauss_paste=" "

90 108
     DO While (LineL.GT.0)

91 109
        READ(IOIN,'(A)') Line

92 110
        Line=AdjustL(Line)

93 111
        LineL=Len_TRIM(Line)

94
        IF (LineL.GT.0) Iat=Iat+1

95
        Idx=Index(Line,'.',BACK=.TRUE.)

96
        Line=ADJUSTL(Line(Idx+1:))

97

98
        IF (LEN_TRIM(Line).GT.0) THEN

99
           Gauss_paste(Iat)=ADJUSTL(TRIM(Line))

112
        IF (LineL.GT.0) THEN

113
           LineUp=Line

114
           Call UpCase(LineUp)

115
           if (LineUp(1:2)=="TV") THEN

116
              FPBC=.TRUE.

117
              IPer=IPer+1

118
              If (Iper>3) THEN

119
                 Call Die("ReadInput Gaussian","Iper>3",Unit=IOOUT)

120
              END IF

121
              NTmp=Index(LineUp," ")

122
              LineUp=LineUp(NTmp:)

123
              Read(LineUp,*) Lat(IPer,1:3)

124
           ELSE

125
              Iat=Iat+1

126
! we search for additional information at the end of the line

127
! for example ONIOM layers

128
! TO detect the end of the line, we use the fact that all reals

129
! should contain a '.', and that we have 3 reals/line.

130
              Idx=Index(Line,'.')

131
              Line=ADJUSTL(Line(Idx+1:))

132
              Idx=Index(Line,'.')

133
              Line=ADJUSTL(Line(Idx+1:))

134
              Idx=Index(Line,'.')

135
              Line=ADJUSTL(Line(Idx+1:))

136
              Idx=Index(Line,' ')

137
              If (Idx>0) THEN

138
                 Line=ADJUSTL(Line(Idx:))

139
                 IF (LEN_TRIM(Line).GT.0) THEN

140
                    Gauss_paste(Iat)=ADJUSTL(TRIM(Line))

141
                 END IF

142
              ELSE

143
                 Gauss_paste(Iat)=""

144
              END IF

145
           END IF

100 146
        END IF

147


101 148
     END DO

102 149

150

103 151
     IF (Iat.NE.Nat) THEN

104
        WRITE(*,*) "I found ", Iat," lines for the geometry instead of ",Nat

105
        WRITE(*,*) "ERROR. STOP"

106
        WRITE(IOOUT,*) "I found ", Iat," lines for the geometry instead of ",Nat

107
        WRITE(IOOUT,*) "ERROR. STOP"

108
        STOP

152
        WRITE(Line,*) "I found ", Iat," lines for the geometry instead of ",Nat

153
        Call Die("ReadInput Gaussian","Line",UNIT=IOOUT)

109 154
     END IF

155

156
  IF (FPBC) THEN

157
     Lat_a(1:3)=Lat(1,1:3)

158
     Lat_b(1:3)=Lat(2,1:3)

159
     Lat_c(1:3)=Lat(3,1:3)

160
     If (IPer>=1) THEN

161
        kaBeg=-1

162
        kaEnd=1

163
     END IF

164
     If (IPer>=2) THEN

165
        kbBeg=-1

166
        kbEnd=1

167
     END IF

168
     If (IPer==3) THEN

169
        kcBeg=-1

170
        kcEnd=1

171
     END IF

172
     If (IPer>3) THEN

173
        Call Die("Readinput_gaussian","Found too many Tv lines !",Unit=IOOUT)

174
     END IF

175
  END IF

176

177

110 178
     ! We now read the last part

111 179
     IF (DEBUG) WRITE(*,*) "Reading Gauss End"

112 180
     !     READ(IOIN,'(A)') Line

......
134 202
        END DO

135 203


136 204
        WRITE(*,*)

137
        WRITE(*,*) 'Comment original'

205
!        WRITE(*,*) '//INFO// Comment original:'

138 206


139 207
        Current => Gauss_comment

140 208
        DO WHILE (ASSOCIATED(Current%next))

......
146 214
        WRITE(*,*) Trim(Gauss_charge)

147 215


148 216
        DO I=1,Nat

149
           WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XyzGeomI(1,1:3,I), TRIM(Gauss_Paste(I))

217
           WRITE(*,'(1X,A10,3(1X,F15.8),1X,A)') Trim(AtName(I)),XyzGeomI(1,1:3,I), TRIM(Gauss_Paste(I))

150 218
        END DO

151 219


152 220
        WRITE(*,*)

3 3
  ! This routines calculates the energy and the gradient of

4 4
  ! a molecule, using Siesta

5 5

6
  use Path_module, only : Nat, renum,Order,OrderInv, Coord, ProgExe

6
  use Path_module, only : Nat, renum,Order,OrderInv, Coord, ProgExe,v_direct,latr

7 7
  use Io_module

8 8

9 9
  IMPLICIT NONE

......
73 73
  logical           :: debug, TChk

74 74
  LOGICAL, SAVE :: first=.true.

75 75

76
  REAL(KREAL) :: Pi

76
  REAL(KREAL), ALLOCATABLE :: X(:),Y(:),Z(:) ! Nat

77
  REAL(KREAL), ALLOCATABLE :: GeomCartLoc(:,:) !Nat,3

77 78

78
  INTEGER(KINT) :: iat, i, n3at,Idx

79
  INTEGER(KINT) :: iat, i, J,Idx

79 80

80 81
  CHARACTER(LCHARS) :: FileIn, FileOut, FileForces,FileE

81 82

......
84 85

85 86
  ! ======================================================================

86 87

87
  Pi=dacos(-1.0d0)

88
  n3at=3*nat

89

90 88
  debug=valid('EGRAD').OR.valid('egrad_siesta')

91 89

92 90
  if (debug) Call Header("Entering Egrad_siesta")

......
103 101
     END IF

104 102
  END IF

105 103

104
  ALLOCATE(GeomCartLoc(Nat,3))

105

106 106
  gradcart=0.

107 107
  E=0.

108 108
  FileIn=Trim(CalcName) // "." // Trim(ISuffix)

......
122 122

123 123
  WRITE(IOTMP,*)

124 124

125
  IF (V_Direct=='DIRECT') THEN

126
     ALLOCATE(X(Nat),Y(Nat),Z(Nat))

127
     X=0.

128
     Y=0.

129
     Z=0.

130
     DO I=1,Nat

131
        DO J=1,3

132
           X(I)=X(I)+GeomCart(I,J)*Latr(1,J)

133
           Y(I)=Y(I)+GeomCart(I,J)*Latr(2,J)

134
           Z(I)=Z(I)+GeomCart(I,J)*Latr(3,J)

135
        END DO

136
     END DO

137
     GeomCartLoc(1:Nat,1)=X(1:Nat)

138
     GeomCartLoc(1:Nat,2)=Y(1:Nat)

139
     GeomCartLoc(1:Nat,3)=Z(1:Nat)

140
     DEALLOCATE(X,Y,Z)

141
  ELSE

142
     GeomCartLoc=GeomCart

143
  END IF

144

125 145
  WRITE(IOTMP,'(1X,A)')  '%block AtomicCoordinatesAndAtomicSpecies'

126 146

127 147
  DO I=1,Nat

......
206 226
     WRITE(*,'(1X,3(1X,F15.8))') GradCart(3*I-2:3*I)

207 227
  END DO

208 228

229
  DeAllocate(GeomCartLoc)

230

209 231
!  Call Die('Egrad_Siesta','Ah ah.',UNIT=IOOUT)

210 232
  if (debug) Call header("Egrad_siesta Over")

211 233

234

235

212 236
  RETURN

213 237

214 238
  ! ======================================================================

src/Path_module.f90 (revision 10)
12 12
  INTEGER(KINT), PARAMETER :: MaxFroz=100

13 13
  REAL(KREAL), PARAMETER :: a0=0.529177249d0

14 14
  REAL(KREAL), PARAMETER :: Unit=1.d0/a0

15
  REAL(KREAL), PARAMETER :: Ang2au=a0, Au2Ang=Unit

15 16
  REAL(KREAL) :: Pi

16 17

17 18
 ! Frozen contains the indices of frozen atoms

......
192 193

193 194
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

194 195
!

196
! Variables for periodic calculations

197
!

198
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

199

200
! Lattice constants

201
  REAL(KREAL) :: lat_a(3), lat_b(3), lat_c(3)

202
! Inverse lattice constants

203
  REAL(KREAL) :: Latr(3,3)

204
! FPBC : True if this is a periodic calculations

205
  LOGICAL :: FPBC

206
! Possible values for ka, kb, kc in VectorPer (and other periodic operations)

207
  INTEGER(KINT) :: kaBeg,kaEnd,kbBeg,kbEnd,kcBeg,kcEnd

208
! Number of periodic directions

209
  INTEGER(KINT) :: IPer

210
! Reference cartesian coordinates

211
  REAL(KREAL), ALLOCATABLE :: XGeomRefPBC(:),YGeomRefPBC(:),ZGeomRefPBC(:) ! Nat

212
! How shall we print the cartesian coordinates ?

213
! if V_Direct='DIRECT' then we use fractional coord (ie divided by unit cell

214
! vectors).

215
  CHARACTER(LCHARS) :: V_direct, V_direct_write

216

217

218
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

219
!

195 220
! Variables for VASP input/output

196 221
!

197 222
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

198 223


199
  REAL(KREAL) :: lat_a(3), lat_b(3), lat_c(3)

200 224
  REAL(KREAL), ALLOCATABLE :: X0_vasp(:),Y0_vasp(:), Z0_vasp(:) ! nat

201
  REAL(KREAL) :: Latr(3,3)

202
  CHARACTER(LCHARS) :: V_direct

203 225
! AutoCart : true if user let PATH determines which atoms should be

204 226
! described in cartesian when COORD=MIXED

205 227
  LOGICAL :: AutoCart

93 93
  DO I=1,Nat

94 94
     If (renum) THEN

95 95
        Iat=Order(I)

96
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),GeomCart(Iat,:), TRIM(Gauss_Paste(I))

96
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8),1X,A)') Trim(AtName(I)),GeomCart(Iat,:), TRIM(Gauss_Paste(I))

97 97
     ELSE

98 98
        Iat=OrderInv(I)

99
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomCart(I,:), TRIM(Gauss_Paste(Iat))

99
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8),1X,A)') Trim(AtName(Iat)),GeomCart(I,:), TRIM(Gauss_Paste(Iat))

100 100
     END IF

101 101
  END DO

102 102

102 102

103 103

104 104
  CHARACTER(132) ::  Line,Line2

105
  INTEGER(KINT) ::  Idx

105
  INTEGER(KINT) ::  Idx,LineL

106 106
  INTEGER(KINT) ::  IAt,I

107 107
  INTEGER(KINT) :: IoRead, ITmp

108
  REAL(KREAL) :: Xtmp, Ytmp, Ztmp

108
  REAL(KREAL) :: B(3),Xtmp, Ytmp, Ztmp

109
  REAL(KREAL) :: Lata,Latb,Latc,Alpha,Beta,Gamma

110
  REAL(KREAL), ALLOCATABLE :: GeomCartTmp(:,:) ! (Nat,3)

111
  REAL(KREAL) :: LatLoc(3,3)

109 112


110 113
  LOGICAL :: Debug

111 114
  LOGICAL :: FSpecies, FCoord

......
312 315
           END IF

313 316
        END IF

314 317

315
        IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESFORMAT",Search,Clean=".-_")) THEN

318
        IF (SearchInput(Siesta_Input,"LATTICE",Search,Clean=".-_")) THEN

319
! This is a pbc calculation

320
           FPBC=.TRUE.

321
           IPer=3

322
           Kabeg=-1

323
           kaEnd=1

324
           kbBeg=-1

325
           kbEnd=1

326
           kcBeg=-1

327
           kcEnd=1

328


329

330
         IF (SearchInput(Siesta_Input,"LATTICECONSTANT",Search,Clean=".-_")) THEN

316 331
           Line=Adjustl(Search%Line)

332
! We discar the label

317 333
           Idx=Index(Line,' ')

318 334
           Line=AdjustL(Line(Idx+1:))

335
! we read the value

336
           Read(Line,*) Siesta_LatticeConstant

337
! We discard the value

338
           Idx=Index(Line,' ')

339
           Line=AdjustL(Line(Idx+1:))

340
           LineL=Len_Trim(Line)

341
! We read the unit

342
           If (LineL>0) THEN

343
              Call UpCase(Line)

344
              SELECT CASE (Line)

345
                CASE ('ANG')

346
                   Siesta_Lat_Unit=1.d0

347
                CASE ('BOHR')

348
                   Siesta_Lat_Unit=a0

349
              END SELECT

350
           ELSE

351
              Siesta_Lat_Unit=1.d0

352
           END IF

353
         ELSE

354
! We have no LatticeConstant but we need one for our calculation

355
            Siesta_LatticeConstant=1.

356
            Siesta_Lat_Unit=1.d0

357
         END IF

358

359
         IF (SearchInput(Siesta_Input,"LATTICEVECTORS",Search,Clean=".-_")) THEN

360
           search => search%next

361
           Line=Adjustl(Search%Line)

362
           Read(Line,*) Lat_a(1:3)

363
           search => search%next

364
           Line=Adjustl(Search%Line)

365
           Read(Line,*) Lat_b(1:3)

366
           search => search%next

367
           Line=Adjustl(Search%Line)

368
           Read(Line,*) Lat_c(1:3)

369

370
         ELSEIF (SearchInput(Siesta_Input,"LATTICEPARAMETERS",Search,Clean=".-_")) THEN

371
           search => search%next

372
           Line=Adjustl(Search%Line)

373
           Read(Line,*) Lata,Latb,Latc,Alpha,Beta,Gamma

374

375
            Lat_a=0.

376
            Lat_b=0.

377
            Lat_c=0.

378
            Lat_a(1)=Lata

379
            Lat_b(1)=LatB*cos(Gamma*Pi/180d0)

380
            Lat_b(2)=LatB*sin(Gamma*Pi/180d0)

381
            Lat_c(1)=Latc*cos(Beta*Pi/180.)

382
            Lat_c(2)=(Latb*Latc*cos(Alpha*Pi/180.)-Lat_b(1)*Lat_c(1))/Lat_b(2)

383
            Lat_c(3)=sqrt(Latc*Latc-Lat_c(1)**2-Lat_c(2)**2)

384
         ELSE

385
! We have a lattice constant but nothing else. We issue a warning just in case,

386
! and we take the defaut values 1. 1. 1. 90. 90. 90.

387
            Line="LatticeConstant found, but not LatticeVectors nor LatticeParameters. Taking 1. 1. 1. 90. 90. 90."

388
            Call Warning("ReadInput_siesta",Line)

389
            Lat_a=(/1.,0.,0./)

390
            Lat_b=(/0.,1.,0./)

391
            Lat_c=(/0.,0.,1./)

392
         END IF

393

394
         IF (debug) THEN

395
            WRITE(*,*) "Lattice vectors before scaling with LatticeConstant"

396
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_a",Lat_a

397
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_b",Lat_b

398
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_c",Lat_c

399
         END IF

400
           Lat_a=Lat_a*Siesta_LatticeConstant*Siesta_Lat_Unit

401
           Lat_b=Lat_b*Siesta_LatticeConstant*Siesta_Lat_Unit

402
           Lat_c=Lat_c*Siesta_LatticeConstant*Siesta_Lat_Unit

403

404
         IF (debug) THEN

405
            WRITE(*,*) "Lattice vectors After scaling with LatticeConstant"

406
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_a",Lat_a

407
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_b",Lat_b

408
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_c",Lat_c

409
         END IF

410

411
      END IF

412

413
         IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESFORMAT",Search,Clean=".-_")) THEN

414
           Line=Adjustl(Search%Line)

415
           Idx=Index(Line,' ')

416
           Line=AdjustL(Line(Idx+1:))

319 417
           Call UpCase(Line)

320 418
           SELECT CASE (Line)

321 419
              CASE ('ANG','NOTSCALEDCARTESIANANG')

322 420
                 Siesta_Unit_Read=1.d0

323 421
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')

422
                 IF (FPBC) THEN

423
                    Siesta_Unit_Read=1.d0

424
                    V_DIRECT="DIRECT"

425
                    Latr(1:3,1)=Lat_a

426
                    Latr(1:3,2)=Lat_b

427
                    Latr(1:3,3)=Lat_c

428
                    B=1.

429
! TO DO: replace by LAPACK

430
                    CALL Gaussj(Latr,3,3,B,1,1)

431
                    LatLoc(1:3,1)=Lat_a

432
                    LatLoc(1:3,2)=Lat_b

433
                    LatLoc(1:3,3)=Lat_c

434
                    Allocate(GeomCartTmp(Nat,3))

435
                    Do I=1,NGeomI

436
                       GeomCartTmp=Reshape(XyzGeomI(I,1:3,1:Nat),(/Nat,3/),ORDER=(/2,1/))

437
! TO DO: shall we replace MatMull by DGEMM ?

438
                       XyzGeomI(I,1:3,1:Nat)=Reshape(MatMul(GeomCartTmp,Latloc),(/3,Nat/),ORDER=(/2,1/))

439
                   END DO

440
                 ELSE

441
                    WRITE(Line2,*) "AtomicCoordinatesFormat=" // Trim(Line) &

442
               //", but there is no information about Lattice parameters"

443
                    Call Die("ReadInput_siesta",Line2)

444
                 END IF

324 445
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')

325 446
                 Siesta_Unit_Read=a0

326 447
                 IF (INPUT=='SIESTA') THEN

327 448
! We have read the coordinates, but not the unit. This is corrected here

328 449
                    XyZGeomI=XyzGeomI*a0

329 450
                 END IF

451
              CASE ('SCALEDCARTESIAN')

452
                 Siesta_Unit_Read=Siesta_LatticeConstant

453
                 If (INPUT=='SIESTA') THEN

454
                    XyzGeomI=XyzGeomI*Siesta_Unit_Read

455
                 END IF

330 456
            END SELECT

457
         ELSE

458
            Siesta_Unit_Read=1.d0

331 459
         END IF

332 460

333 461
        IF (SearchInput(Siesta_Input,"ATOMICCOORFORMATOUT",Search,Clean=".-_")) THEN

......
339 467
              CASE ('ANG','NOTSCALEDCARTESIANANG')

340 468
                 Siesta_Unit_Write=1.d0

341 469
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')

470
                 Siesta_Unit_Write=1.d0

471
                 V_DIRECT_Write='DIRECT'

342 472
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')

343
                 Siesta_Unit_Write=Unit

473
                 Siesta_Unit_Write=Unit ! 1/a0

344 474
            END SELECT

345 475
         END IF

346

347
         IF (SearchInput(Siesta_Input,"LATTICECONSTANT",Search,Clean=".-_")) THEN

348
           Line=Adjustl(Search%Line)

349
! We discar the label

350
           Idx=Index(Line,' ')

351
           Line=AdjustL(Line(Idx+1:))

352
! we read the value

353
           Read(Line,*) Siesta_LatticeConstant

354
! We discard the value

355
           Idx=Index(Line,' ')

356
           Line=AdjustL(Line(Idx+1:))

357
! We read the unit

358
           Call UpCase(Line)

359
           SELECT CASE (Line)

360
              CASE ('ANG')

361
                 Siesta_Lat_Unit=1.d0

362
              CASE ('BOHR')

363
                 Siesta_Lat_Unit=Unit

364
            END SELECT

365
! for now!

366
            Call Die('ReadInput_siesta:LatticeConstant',"For now, periodic calculations are NOT possible in Opt'n Path")

367

368
         END IF

369

370
         IF (SearchInput(Siesta_Input,"LATTICEVECTORS",Search,Clean=".-_")) THEN

371

372
            Call Die('ReadInput_siesta:LatticeVectors',"For now, periodic calculations are NOT possible in Opt'n Path")

373

374
         END IF

375

376
         IF (SearchInput(Siesta_Input,"LATTICEPARAMETERS",Search,Clean=".-_")) THEN

377

378
            Call Die('ReadInput_siesta:LatticeParameters',"For now, periodic calculations are NOT possible in Opt'n Path")

379

380
         END IF

381

382 476
         IF (SearchInput(Siesta_Input,"SUPERCELL",Search,Clean=".-_")) THEN

383

384
            Call Die('ReadInput_siesta:SuperCell',"SuperCell  NOT possible in Opt'n Path")

385

477
            Line=Adjustl(search%line)

478
            Call Upcase(Line)

479
            Call CleanString(Line,".-_")

480
            If (Index(Line,"%BLOCK")/=0) THEN

481
               Call Die('ReadInput_siesta:SuperCell',"SuperCell  NOT possible in Opt'n Path")

482
            END IF

386 483
         END IF

387 484

388 485
!         Call Die('Readinput_siesta:fin',' Lecture finie')

8 8
       prog,NCart,XyzGeomF, XyzGeomI, renum, &

9 9
       GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &

10 10
       , BprimT,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &

11
      , OrderInv, XPrimitiveF

11
      ,Order, OrderInv, XPrimitiveF,FPBC,XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC

12 12

13 13
  ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)

14 14
  ! allocated in Path.f90

......
43 43
  REAL(KREAL), ALLOCATABLE :: GradCart(:)

44 44
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)

45 45
  REAL(KREAL) ::  Pi

46
  REAL(KREAL) :: MRot(3,3), Rmsd

46 47

47 48
  INTEGER(KINT) :: iat, i, j, IBeg,Idx

48 49

......
92 93
  CASE ('ZMAT')

93 94
     ! In order to avoid problem with numbering and linear angles/diehedral, we convert the

94 95
     ! zmat into cartesian coordinates

95
     ! A remplacer par Int2Cart :-)

96 96
     Call Int2Cart(nat,indzmat,Geom,GeomCart)

97
! In case of PBC calculations, we have to re-orient the geometry onto the user initial geometry

98
     If (FPBC) THEN

99
! we align this geometry on the initial one

100
        if (debug) THEN

101
           WRITe(*,*) "We are orientating..."

102
           WRITE(*,*) "Geom before orientation"

103
           WRITE(*,*) Nat

104
           WRITE(*,*) ""

105
           DO I=1,Nat

106
              If (renum) Iat=Order(I)

107
              WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomCart(I,1),GeomCart(I,2),GeomCart(I,3)

108
           END DO

109
        END IF

110

111
        IF (Renum) THEN

112
           DO I=1,Nat

113
              Iat=Order(I)

114
              X(I)=GeomCart(Iat,1)

115
              Y(I)=GeomCart(Iat,2)

116
              Z(I)=GeomCart(Iat,3)

117
           END DO

118
        ELSE

119
           X(1:Nat)=GeomCart(1:Nat,1)

120
           Y(1:Nat)=GeomCart(1:Nat,2)

121
           Z(1:Nat)=GeomCart(1:Nat,3)

122
        END IF

123

124
           WRITE(*,*) "Geom before orientation after reorderring"

125
           WRITE(*,*) Nat

126
           WRITE(*,*) ""

127
           DO I=1,Nat

128
!              Iat=I

129
!              If (Renum) Iat=Order(I)

130
              WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),X(I),Y(I),Z(I)

131
           END DO

132
           WRITE(*,*) "Ref Geom"

133
           WRITE(*,*) Nat

134
           WRITE(*,*) ""

135
           DO I=1,Nat

136
              WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XGeomRefPBC(I),YGeomRefPBC(I),ZGeomRefPBC(I)

137
           END DO

138
        END IF

139

140
        Call CalcRmsd(Nat,XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC, &

141
             X,Y,Z, MRot, Rmsd,.TRUE.,.TRUE.)

142

143
        If (DEBUG) THEN

144
           WRITE(*,*) "Geom AFTER orientation"

145
           WRITE(*,*) Nat

146
           WRITE(*,*) ""

147
           DO I=1,Nat

148
              WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),X(I),Y(I),Z(I)

149
           END DO

150
        END IF

151


152
        If (Renum)  THEN

153
           Do I=1,Nat

154
              Iat=orderInv(I)

155
              GeomCart(Iat,1)=X(I)

156
              GeomCart(Iat,2)=Y(I)

157
              GeomCart(Iat,3)=Z(I)

158
           END DO

159
        END IF

97 160
  CASE ('BAKER')

98 161
     XPrimRef=XPrimitiveF(IGeom,:)

99 162
     IF (Flag_Opt_Geom) THEN

......
184 247
  !  WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3)

185 248
  !END DO

186 249

250

251

187 252
  SELECT CASE (Prog)

188 253
  CASE ('GAUSSIAN')

189 254
! we call the Gaussian routine.

33 33

34 34
  if (debug) Call Header("Entering ReadInput_Vasp")

35 35

36
  FPBC=.TRUE.

37

36 38
 if (Input/="VASP") THEN

37 39

38 40
     ! Input was not Vasp, so many parameters are missing like lattice

......
138 140
        Latr(1:3,2)=Lat_b

139 141
        Latr(1:3,3)=Lat_c

140 142
        B=1.

143
! TO DO: replace by LAPACK

141 144
        CALL Gaussj(Latr,3,3,B,1,1)

142 145
     ELSE

143 146
        Latr=0.

......
180 183
        Vasp_Title=Trim(Line) // " " // Trim(adjustL(Vasp_Title))

181 184
        if (debug) WRITE(*,*) "DBG MAIN, L809, VMD=T, Vasp_Title=",TRIM(Vasp_Title)

182 185
     END IF

183
     Call CheckPeriodicBound

184 186

185 187
 if (debug) Call Header("Exiting ReadInput_Vasp")

186 188

src/Path.f90 (revision 10)
280 280

281 281
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

282 282

283
  LOGICAL :: Debug, Fini,Flag_Opt_Geom

283
  LOGICAL :: Debug, Fini,Flag_Opt_Geom,TChk

284 284

285 285
!  INTEGER(KINT), EXTERNAL :: Iargc

286 286
  INTEGER(KINT), EXTERNAL ::  ConvertNumAt

......
831 831
     If (ProgExe=="EMPTY") ProgExe="VASP"

832 832
     UnitProg="eV"

833 833
     ConvE=eV2kcal

834
     FPBC=.TRUE.

835
     IPer=3

836
     Kabeg=-1

837
     kaEnd=1

838
     kbBeg=-1

839
     kbEnd=1

840
     kcBeg=-1

841
     kcEnd=1

834 842
    CASE ("TM","TURBOMOLE")

835 843
     Prog="TURBOMOLE"

836 844
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"

837 845
     UnitProg="au"

838 846
     ConvE=au2kcal

847
     FPBC=.FALSE.

839 848
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")

840 849
     ProgExe=""

841 850
     UnitProg="au"

......
1192 1201

1193 1202
  Call ReadInput

1194 1203

1204
  if (FPBC)  THEN

1205
     Allocate(XGeomRefPBC(Nat),YGeomRefPBC(Nat),ZGeomRefPBC(Nat))

1206
     XGeomRefPBC(1:Nat)=XyzGeomI(1,1,1:Nat)

1207
     YGeomRefPBC(1:Nat)=XyzGeomI(1,2,1:Nat)

1208
     ZGeomRefPBC(1:Nat)=XyzGeomI(1,3,1:Nat)

1209

1210
     Call CheckPeriodicBound

1211
  END IF

1212

1195 1213
  IF (COORD=="MIXED") Call TestCart(AutoCart)

1196 1214

1197 1215
  SELECT CASE(NFroz)

......
2345 2363
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)

2346 2364
  END IF

2347 2365

2366
  If (FPBC)   DeAllocate(XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC)

2348 2367
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt

2349 2368

2350 2369
  Close(IOIN)

src/NoComment.f90 (revision 10)
1
      SUBROUTINE noComment(String)

2

3
! This subroutine suppress all comments from a string

4
! comments can start with ! or #

5

6
        IMPLICIT NONE

7

8
        integer, parameter :: KINT = kind(1)

9

10
        CHARACTER(*), INTENT(INOUT) :: String

11
        INTEGER(KINT) :: i

12

13
        I=Scan(String,"!#")

14
        IF (I/=0) THEN

15
           String=String(:I-1)

16
        END IF

17

18
      END SUBROUTINE NoComment

19

src/NoString.f90 (revision 10)
1
      SUBROUTINE noString(String)

2

3
! This subroutine suppress the content of string variable

4
! strings can be delimited by ' or "

5
! example:

6
! Progexe='./siesta', -> Progexe=,

7

8
        use VarTypes

9

10
        IMPLICIT NONE

11

12
        CHARACTER(*), INTENT(INOUT) :: String

13
        CHARACTER(VLCHARS) :: Result

14
        INTEGER(KINT) :: i

15
        CHARACTER(2), PARAMETER :: Del='"' // "'"

16

17
        DO WHILE (scan(String,Del)/=0)

18
           I=Scan(String,Del)

19
           Result=String(:I-1)

20
           String=String(I+1:)

21
           I=Scan(String,Del)

22
           String=TRIM(Result) // String(I+1:)

23
        END DO

24

25
      END SUBROUTINE NoString

26

src/CalcCnct.f90 (revision 10)
30 30

31 31
SUBROUTINE CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)

32 32

33
  use Path_module, only : NMaxL, max_Z,Nom,Prog, KINT, KREAL

33
  use Path_module, only : NMaxL, max_Z,Nom,Prog, KINT, KREAL, FPBC, &

34
       kaBeg,kaEnd,kbBeg,kbEnd,kcBeg,kcEnd

34 35


35 36
  IMPLICIT NONE

36 37

......
64 65
  END DO

65 66

66 67
  if (debug) THEN

67
     WRITE(*,*) 'CalcCnct'

68
     WRITE(*,*) 'CalcCnct : covalent radii used'

68 69
     DO iat=1,na

69 70
        i=atome(iat)

70
        WRITE(*,*) Nom(I),I,r_cov(i),r_cov(i)*fact

71
        WRITE(*,*) Nom(I),I,r_cov(i)*fact

71 72
     END DO

73
     WRITE(*,*) 'Coordinates'

74
     DO iat=1,na

75
        i=atome(iat)

76
        WRITE(*,*) Nom(I),x(iat),y(iat),z(iat)

77
     END DO

72 78
  END IF

73 79

74
  IF (PROG/="VASP") THEN

80
  IF (FPBC) THEN

75 81
     DO i=1,na

76 82
        NbLi=LIAISONS(i,0)

77 83
        DO j=i+1,na

78
           CALL vecteur(j,i,x,y,z,vx,vy,vz,dist)

79
           dist=dist/fact

80
           distth=(r_cov(atome(i))+r_cov(atome(j)))/100.

81
!           if (debug) WRITE(*,*) atome(i),atome(j),dist,distth

82
           if (dist.le.distth) THEN

84
           Bound=.FALSE.

85
           DO ka=kaBeg,kaEnd

86
              DO Kb=kbBeg,kbEnd

87
                 DO Kc=kcBeg,kcEnd

88
                    CALL VectorPer(j,i,ka,kb,kc,x,y,z,vx,vy,vz,dist)

89
                    dist=dist/fact

90
                    distth=(r_cov(atome(i))+r_cov(atome(j)))/100.

91
!                    if (debug) WRITE(*,*) atome(i),atome(j),dist,distth

92
                    if (dist.le.distth) Bound=.TRUE.

93
                 END DO

94
              END DO

95
           END DO

96
           IF (Bound) THEN

97
              if (debug) WRITE(*,*) "Adding a bond between:",i,j

83 98
              NbLi=NbLi+1

84 99
              LIAISONS(i,NbLi)=j;

85 100
              NBlj=LIAISONS(j,0)+1

......
93 108
     DO i=1,na

94 109
        NbLi=LIAISONS(i,0)

95 110
        DO j=i+1,na

96
           Bound=.FALSE.

97
           DO ka=-1,1

98
              DO Kb=-1,1

99
                 DO Kc=-1,1

100
                    CALL VectorPer(j,i,ka,kb,kc,x,y,z,vx,vy,vz,dist)

101
                    dist=dist/fact

102
                    distth=(r_cov(atome(i))+r_cov(atome(j)))/100.

103
!                    if (debug) WRITE(*,*) atome(i),atome(j),dist,distth

104
                    if (dist.le.distth) Bound=.TRUE.

105
                 END DO

106
              END DO

107
           END DO

108
           IF (Bound) THEN

109
              if (debug) WRITE(*,*) "Adding a bond between:",i,j

111
           CALL vecteur(j,i,x,y,z,vx,vy,vz,dist)

112
           dist=dist/fact

113
           distth=(r_cov(atome(i))+r_cov(atome(j)))/100.

114
!           if (debug) WRITE(*,*) atome(i),atome(j),dist,distth

115
           if (dist.le.distth) THEN

110 116
              NbLi=NbLi+1

111 117
              LIAISONS(i,NbLi)=j;

112 118
              NBlj=LIAISONS(j,0)+1

src/CheckPeriodicBound.f90 (revision 10)
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.

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 6

7 7

8 8
  Use Path_module

......
15 15
  REAL(KREAL), ALLOCATABLE :: GeomRef(:,:), GeomTmp(:,:) ! (3,Nat)

16 16

17 17
  REAL(KREAL) :: Latrloc(3,3),Bloc(3),Xt,Yt,Zt

18
  REAL(KREAL) :: NormA, NormB,Xrel,Yrel,ZRel,NormC

19
  REAL(KREAL) :: Prod,LatrA(3),LatrB(3)

20
  REAL(KREAL) :: V(3)

21

18 22
  LOGICAL :: Debug

19 23

20 24
  INTERFACE

......
22 26
       CHARACTER(*), intent(in) :: string

23 27
       logical                  :: isValid

24 28
     END function VALID

29

30

31
    SUBROUTINE die(routine, msg, file, line, unit)

32

33
      Use VarTypes

34
      Use io_module

35

36
      implicit none

37
      character(len=*), intent(in)           :: routine, msg

38
      character(len=*), intent(in), optional :: file

39
      integer(KINT), intent(in), optional      :: line, unit

40

41
    END SUBROUTINE die

42

25 43
  END INTERFACE

26 44

27 45
  debug=valid('checkperiodicbound')

28
if (debug) WRITE(*,*) "========================= Entering CheckPeriodicBound ==========="

46
  if (debug) call header("Entering CheckPeriodicBound")

29 47

30 48
  ALLOCATE(GeomRef(3,Nat), GeomTmp(3,Nat))

31 49

32
  V_direct=adjustl(Vasp_direct)

33
  Call upcase(V_direct)

50
  ! The test is easy in DIRECT space, so we will convert the geometries

51
  ! into this space.

34 52

35
! The test is easy in DIRECT space, so we will convert the geometries

36
! into this space.

37
! We might optimize this as Read_geom converts geometries into NORMAL

38
! coordinates... so it might be that we are converting them back to direct here...

53
  SELECT CASE (Iper)

54
  CASE(1)

55
     NormA=DOT_PRODUCT(Lat_a,Lat_a)

56
     DO IGeom=2,NGeomI

57
        DO I=1,NAt

58
           DO J=1,3

59
              v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)

60
           END DO

61
           XRel=DOT_PRODUCT(Lat_a,v)/NormA

62
           if (abs(Xrel).GE.BoxTol) THEN

63
              If (debug) THEN

64
                 WRITE(*,'("Atom ",I5," moved because XRel=",F10.6)') I,XRel

65
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)

66
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)

67
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Lat_a*Sign(1.0d0,XRel)

68
              END IF

69
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Lat_a*Sign(1.0d0,XRel)

70
           END IF

71
        END DO

72
     END DO

73
  CASE (2)

74
     NormA=DOT_PRODUCT(Lat_a,Lat_a)

75
     NormB=DOT_PRODUCT(Lat_b,Lat_b)

76
     Prod=DOT_PRODUCT(Lat_a,Lat_b)

77
     LatrA=(NormB*Lat_a-Prod*Lat_b)/(NormA*NormB-Prod**2)

78
     LatrB=(NormA*Lat_b-Prod*Lat_a)/(NormA*NormB-Prod**2)

39 79

80
     DO IGeom=2,NGeomI

81
        DO I=1,NAt

82
           DO J=1,3

83
              v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)

84
           END DO

85
           XRel=DOT_PRODUCT(Latra,v)

86
           YRel=DOT_PRODUCT(Latrb,v)

87
           if (abs(Xrel).GE.BoxTol) THEN

88
              If (debug) THEN

89
                 WRITE(*,'("Atom ",I5," moved because XRel=",F10.6)') I,XRel

90
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)

91
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)

92
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Latra*Sign(1.0d0,XRel)

93
              END IF

94
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Latra*Sign(1.0d0,XRel)

95
           END IF

96
           if (abs(Yrel).GE.BoxTol) THEN

97
              If (debug) THEN

98
                 WRITE(*,'("Atom ",I5," moved because YRel=",F10.6)') I,YRel

99
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)

100
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)

101
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Latrb*Sign(1.0d0,YRel)

102
              END IF

103
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Latrb*Sign(1.0d0,YRel)

104
           END IF

105
        END DO

106
     END DO

107

108
  CASE(3)

109

110

111

112
!      if (debug) THEN

113
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_a:",Lat_a

114
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_b:",Lat_b

115
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_c:",Lat_c

116
!         WRITE(*,'(1X,A,3(1X,F10.6))') "BocTol:",BoxTol

117

118
!      NormA=DOT_PRODUCT(Lat_a,Lat_a)

119
!      NormB=DOT_PRODUCT(Lat_b,Lat_b)

120
!      Prod=DOT_PRODUCT(Lat_a,Lat_b)

121
!      LatrA=(NormB*Lat_a-Prod*Lat_b)/(NormA*NormB-Prod**2)

122
!      LatrB=(NormA*Lat_b-Prod*Lat_a)/(NormA*NormB-Prod**2)

123
! !     NormA=DOT_PRODUCT(Latra,Latra)

124
! !     NormB=DOT_PRODUCT(Latrb,Latrb)

125

126
!         DO IGeom=2,NGeomI

127
!            DO I=1,NAt

128
!               DO J=1,3

129
!                  v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)

130
!               END DO

131
!               XRel=DOT_PRODUCT(Latra,v)

132
!               YRel=DOT_PRODUCT(Latrb,v)

133
!               If (abs(XRel).GT.BoxTol)  WRITE(*,*) "Xrel too big:",IGeom,I,Xrel

134
!               If (abs(YRel).GT.BoxTol)  WRITE(*,*) "Yrel too big:",IGeom,I,Yrel

135
! !               WRITE(*,*) "Xrel :",IGeom,I,Xrel

136
! !               WRITE(*,*) "Yrel :",IGeom,I,Yrel

137
! !               WRITE(*,*) "Zrel :",IGeom,I,Zrel

138

139
!            END DO

140
!         END DO

141
!      END IF

142

40 143
     Latrloc(1:3,1)=Lat_a

41 144
     Latrloc(1:3,2)=Lat_b

42 145
     Latrloc(1:3,3)=Lat_c

43 146
     Bloc=1.

147
     ! TODO: Replace by LAPACK subroutine

44 148
     CALL Gaussj(Latrloc,3,3,Bloc,1,1)

45 149

46
 ! We put first geometry as the reference

150
     ! We put first geometry as the reference

47 151
     GeomRef=0.

48 152
     DO I=1,Nat

49 153
        DO K=1,3

......
53 157
        END DO

54 158
     END DO

55 159

56
  DO IGeom=2,NGeomI

57
     GeomTmp=0.

58
     DO I=1,Nat

59
        DO K=1,3

60
           DO J=1,3

61
              GeomTmp(K,I)=GeomTmp(K,I)+XyzGeomI(IGeom,J,I)*LatrLoc(K,J)

160
     DO IGeom=2,NGeomI

161
        GeomTmp=0.

162
        DO I=1,Nat

163
           DO K=1,3

164
              DO J=1,3

165
                 GeomTmp(K,I)=GeomTmp(K,I)+XyzGeomI(IGeom,J,I)*LatrLoc(K,J)

166
              END DO

167
              if (abs(GeomTmp(K,I)-GeomRef(K,I)).GE.BoxTol) THEN

168
                 If (debug) &

169
                      WRITE(*,'("Coord ",I1," has been changed for atom ",I5," by ",F4.0," ref:",F9.6," old:",F9.6," new:",F9.6)') &

170
                      K,I,Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)),&

171
                      GeomRef(K,I),GeomTmp(K,I), &

172
                      GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I))

173

174
                 GeomTmp(K,I)=GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I))

175
              END IF

62 176
           END DO

63
           if (abs(GeomTmp(K,I)-GeomRef(K,I)).GE.BoxTol) THEN

64
              If (debug) &

65
                  WRITE(*,'("Coord ",I1," has been changed for atom ",I5," by ",F4.0," ref:",F9.6," old:",F9.6," new:",F9.6)') &

66
                   K,I,Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)),&

67
                   GeomRef(K,I),GeomTmp(K,I), &

68
                   GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I))

69
              GeomTmp(K,I)=GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I))

70
           END IF

177
           Xt=GeomTmp(1,I)

178
           Yt=GeomTmp(2,I)

179
           Zt=GeomTmp(3,I)

180
           XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1)

181
           XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2)

182
           XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3)

71 183
        END DO

72
        Xt=GeomTmp(1,I)

73
        Yt=GeomTmp(2,I)

74
        Zt=GeomTmp(3,I)

75
        XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1)

76
        XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2)

77
        XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3)

78 184
     END DO

79
  END DO

80 185

81
if (debug) WRITE(*,*) "========================= Exiting CheckPeriodicBound ==========="

186
  CASE DEFAULT

187
     Call Die("CheckPeriodicBound"," Iper unknown")

188
  END SELECT

189

190

191
  if (debug) WRITE(*,*) "========================= Exiting CheckPeriodicBound ==========="

82 192
END SUBROUTINE CheckPeriodicBound

83 193

84 194

Also available in: Unified diff