Revision 10

src/egrad_mopac.f90 (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

  
src/ReadInput_mopac.f90 (revision 10)
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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
src/ReadInput_gaussian.f90 (revision 10)
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(*,*) 
src/egrad_siesta.f90 (revision 10)
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
src/egrad_gaussian.f90 (revision 10)
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

  
src/ReadInput_siesta.f90 (revision 10)
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')
src/Egrad.f90 (revision 10)
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.
src/ReadInput_vasp.f90 (revision 10)
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