Statistiques
| Révision :

root / utils / Xyz2Path.f90 @ 8

Historique | Voir | Annoter | Télécharger (29 ko)

1
          program Xyz2Path
2
! This programs reads a XYZ file and converts it into distances,
3
! valence angle and dihedral angles.
4
! It prints them as a function of the irc distance...
5
!-----------------------------------------------
6
! Input: name of the XYZ File
7
! it also needs a file call list which has the following structure:
8
! one line contains the type of the value you want to follow, it can be
9
! b  for a Bond distance
10
! a for an angle
11
! d for a dihedral
12
! this descriptor is followed by the number of the atoms involved !
13
! a typical file can be:
14
! b  1  2
15
! b 2  3
16
! a 1 2 3
17
!----------------------------------------------
18
! Ouput: A files call Scan.dat
19
! wich contains in the first lines the input file (as a reminder)
20
! and then for each step the wanted values
21
!------------------------------------------------
22
! Second version also reads the energy (as to be written after E= on
23
! the comment line)
24
!------------------------------------------------
25
! Third version contains a new command: c for Center of Mass
26
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
27
! v 3.1 the c command now creates the center of mass... and allows
28
! people to do whatever they want with it...
29
! Syntax: c NbAt ListAt
30

    
31

    
32
          IMPLICIT NONE
33

    
34
          INCLUDE "Xyz2Path.param"
35

    
36
          character(40) :: f1
37
          REAL(8) ::   geos(3,maxnat), geos1(3,maxnat)
38
          character(33) :: fmt
39
          character(3) :: atoms(maxnat)
40
          character(5) :: Type
41
          Character(120) :: line
42
          INTEGER(4) :: NbPrint
43
          REAL(8) ::   AU2PS,Pi
44
          REAL(8) ::  Mass(MaxNat), Ener, Conv, Ds, s
45
          INTEGER(4) :: At1,At2,At3,At4,IOOUT,Iat
46
          INTEGER(4) ::  IArg, I, NNN, Ng, J
47

    
48
          INTEGER(4) :: Nat,NbDist, NbAngle, NbDie,NbCOM
49
          INTEGER(4) :: At1B(MaxNat),At2B(MaxNat)
50
          INTEGER(4) :: At1A(MaxNat),At2A(MaxNat),At3A(MaxNat)
51
          INTEGER(4) :: At1D(MaxNat),At2D(MaxNat),At3D(MaxNat),At4D(MaxNat)
52
          INTEGER(4) :: AtCom(0:MaxNat,MaxNat)
53
          REAL(8) ::    VB(MaxNat),VA(MaxNat),VD(MaxNat),VCOM(MaxNat)
54

    
55
          REAL(8) :: MRot(3,3), Rmsd
56
          LOGICAL FExist,FRot,FAlign,Debug,DebugRMSD
57

    
58
          INTEGER(4) :: ConvertNumAt
59
          external ConvertNumAt
60

    
61
          COMMON /Indices/Nat,NbDist,NbAngle,NbDie,NbCom, At1B,At2B, &
62
              At1A,At2A,At3A,At1D,At2D,At3D,At4D,AtCom
63
          COMMON /Values/VB,VA,VD,VCom
64
          COMMON /Const/AU2ps,Pi
65
 
66

    
67

    
68
  REAL(8) :: MassAt(0:86)=(/0.0D0,1.0078D0,          4.0026D0,      &
69
                  7.0160D0, 9.0122D0,11.0093D0,                          &
70
                 12.0000D0,14.0031D0,15.9949D0,18.9984D0,19.9924D0,      &
71
                 22.9898D0,23.9850D0,26.9815D0,                          &
72
                 27.9769D0,30.9738D0,31.9721D0,34.9688D0,39.9624D0,      &
73
                 39.0983D0,40.08D0,                                      &
74
                   44.9559D0, 47.88D0, 50.9415D0, 51.996D0, 54.9380D0,   &
75
                   55.847D0, 58.9332D0, 58.69D0, 63.546D0, 65.39D0,      &
76
                 69.72D0,72.59D0,74.9216D0,78.96D0,79.904D0,83.80D0,     &
77
                 85.4678D0,87.62D0,88.9059D0,91.224D0,92.9064D0,         &
78
         95.94D0,98D0,101.07D0,102.906D0,106.42D0,107.868D0,112.41D0,    &
79
         114.82D0,118.71D0,121.75D0,127.60D0,126.905D0,131.29D0,         &
80
!     6           'CS','BA',
81
                 132.905D0,137.34D0,                                      &
82
!     6       'LA',
83
!               'CE','PR','ND','PM','SM','EU','GD',
84
!               'TB','DY','HO', 'ER','TM','YB','LU',   
85
            138.91D0,                                                      &
86
            140.12D0, 130.91D0, 144.24D0,147.D0,150.35D0, 151.96D0,157.25D0,  &
87
          158.924D0, 162.50D0, 164.93D0, 167.26D0,168.93D0,173.04D0,174.97D0, &
88
!     6                'HF','TA',' W','RE','OS','IR','PT',
89
!                      'AU','HG',
90
!     6                                 'TL','PB','BI','PO','AT','RN'/
91
         178.49D0, 180.95D0, 183.85D0, 186.2D0, 190.2D0, 192.2D0, 195.09D0,  &
92
             196.97D0, 200.59D0,                                             &
93
             204.37D0, 207.19D0,208.98D0,210.D0,210.D0,222.D0 /)
94

    
95
          Debug=.FALSE.
96
          DebugRMSD=.FALSE.
97

    
98
          AU2PS=1./41341.37
99
          NbPrint=100
100
          Pi=dacos(-1.d0)
101
          IOOUT=13
102
          Conv=1.
103
          FRot=.TRUE.
104
          FAlign=.TRUE.
105
          NbDist=0
106
          NbAngle=0
107
          NbDie=0
108

    
109
          iarg=command_argument_count()
110
          if (iarg.lt.1) then
111
                 write(*,*) 'XYZ filename:'
112
              read(*,'(a)') f1
113
            else
114
              call getarg(1,f1)
115
          endif
116

    
117
          open(13,file='Scan.dat')
118

    
119
          INQUIRE(File='list',Exist=FExist)
120
          if (.NOT.FExist) THEN
121
             WRITE(*,*) "No file 'list', just printing Energy"
122
          END IF
123

    
124
          open(11,file=f1)
125
          call rtraitem(11,nnn,ener,geos1,atoms)
126
          close(11)
127

    
128
         DO I=1,nnn
129
            Iat=ConvertNumAt(atoms(I))
130
            Mass(I)=MassAt(Iat)
131
!           write(*,*) I,Atoms(I),Mass(I),geos1(:,I)
132
         END DO
133

    
134
         Nat=nnn
135

    
136
         if (FExist) THEN
137
         open(14,file='list')
138
         Type="d"
139
         DO WHILE (Type.ne."E")
140
          CALL READLINE(14,Type,Line)
141
!          WRITE(*,*) Line,Type
142
          if (Type.eq."b") THEN
143
             NbDist=NbDist+1
144
             READ(Line,*) At1,At2
145
             At1B(NbDist)=At1
146
             At2B(NbDist)=At2
147
             WRITE(13,'("# b ",2(A3,I3))') Atoms(At1),At1,Atoms(At2),At2
148
             WRITE(*,'("# b ",2(A3,I3))') Atoms(At1),At1,Atoms(At2),At2
149
          END IF
150
          if (Type.eq."a") THEN
151
             NbAngle=NbAngle+1
152
             READ(Line,*) At1,At2,At3
153
             At1A(NbAngle)=At1
154
             At2A(NbAngle)=At2
155
             At3A(NbAngle)=At3
156
            WRITE(13,'("# a ",3(A3,I3))') Atoms(At1),At1,Atoms(At2), &
157
                 At2, Atoms(At3),At3
158
            WRITE(*,'("# a ",3(A3,I3))') Atoms(At1),At1,Atoms(At2), &
159
                At2, Atoms(At3),At3
160
          END IF
161
          if (Type.eq."d") THEN
162
             NbDie=NbDie+1
163
             READ(Line,*) At1,At2,At3,At4
164
             At1D(NbDie)=At1
165
             At2D(NbDie)=At2
166
             At3D(NbDie)=At3
167
             At4D(NbDie)=At4
168
            WRITE(13,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2,  &
169
               Atoms(At3),At3,Atoms(At4),At4
170
            WRITE(*,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2,   &
171
               Atoms(At3),At3,Atoms(At4),At4
172

    
173
          END IF
174
          if (Type.eq."c") THEN
175
            NbCOM=NbCOM+1
176
            READ(Line,*) At1,(AtCom(j,NbCOM),j=1,At1)
177
            AtCom(0,NbCOM)=At1
178
            Atoms(nat+NbCom)="G"
179
            WRITE(13,'("# c ",I4,20(A3,I3))') At1,                  &
180
               (Atoms(AtCom(i,NbCoM)),AtCom(i,NbCOM),i=1,At1)
181
            WRITE(*,'("# c ",I4,20(A3,I3))') At1,                   &
182
               (Atoms(AtCom(i,NbCoM)),AtCom(i,NbCOM),i=1,At1)
183
          END IF
184

    
185
         END DO
186
         END IF
187

    
188

    
189

    
190
          fmt='(   (1X,F12.5),1X,F15.6)'
191
          write(fmt(2:4),'(i3)') NbDist+NbAngle+NbDie+1
192
!           write(*,*) nat3
193
!          write(*,*) 'fmt:',fmt
194

    
195
          ng=1
196

    
197
         s=0.
198

    
199
          open(11,file=f1)
200

    
201
10           call rtraitem(11,nnn,ener,geos,atoms)
202
!           WRITE(*,*) nnn
203
           if (nnn.gt.0) then
204

    
205
          if (debug) THEN
206
             WRITE(*,*) "Debug ref geom: geos1"
207
             DO I=1,nnn
208
                WRITE(*,'(1X,I5,1X,F5.1,3(1X,F12.6))') I,mass(I),geos1(:,I)
209
             END DO
210
          END IF
211

    
212
          call CalcRmsd(nnn, geos1, geos,MRot,rmsd,FRot,FAlign,debugRMSD)
213

    
214
          if (debug) THEN
215
             WRITE(*,*) "Debug Geos aligned on geos1,rmsd=",rmsd
216
             DO I=1,nnn
217
                WRITE(*,'(1X,I5,1X,F5.1,3(1X,F12.6))') I,mass(I),geos(:,I)
218
             END DO
219
          END IF
220

    
221
         ds=0.
222
           DO I=1,nnn
223
             DO J=1,3
224
             ds=ds+mass(I)*(geos(J,I)-geos1(J,I))**2
225
             geos1(J,I)=Geos(J,I)
226
              END DO
227
             if (debug) write(*,*) I,ds
228
            END DO
229
            if (debug) WRITE(*,*) "Debug, ds",sqrt(ds)
230
           s=s+sqrt(ds)
231

    
232
! We convert coordinates in a0 into Angstroem
233
!              write(*,*) "Analyse !"
234
            if (FExist) THEN
235
            Call Analyse(geos)
236

    
237
             write(IOOUT,fmt) s,(VB(j)/Conv,j=1,NbDist), &
238
                 (VA(j),j=1,NbAngle), &
239
                 (VD(j),j=1,NbDie),ener
240
             write(*,fmt) s,(VB(j)/Conv,j=1,NbDist), &
241
                 (VA(j),j=1,NbAngle), &
242
                 (VD(j),j=1,NbDie),ener
243
             ELSE
244
                        write(IOOUT,fmt) s,ener
245
             write(*,fmt) s,ener
246
            END IF 
247
           ng=ng+1
248
            goto 10
249
           endif
250
           WRITE(*,*) 'Found ',ng-1,' geometries'
251
           close(11)
252
         end
253

    
254
!--------------------------------------------------------------
255
          subroutine rtraitem(ifil,nnn,E,tab,atoms)
256
!           implicit REAL(8) :: (a-h,o-z)
257
            IMPLICIT NONE
258

    
259
           character(120) :: line
260
           character(3) :: Atoms(*)
261
           INTEGER(4) :: nnn, IFil, Idx, I, J
262
           REAL(8) ::  Tab(3,*), E
263

    
264

    
265
           read(ifil,*,err=99,end=99) nnn
266
           read(ifil,'(A)') line
267
           idx=index(line,'E')
268
!           WRITE(*,*) 'idx,line',idx,line
269
           if (idx/=0) THEN
270
             line=line(idx+2:)
271
              idx=index(line,"=")
272
              if (idx/=0) line=line(idx+1:)
273
              idx=index(line,":")
274
              if (idx/=0) line=line(idx+1:)
275
              read(line,*) E
276
           ELSE
277
             E=0.
278
           END IF
279
           
280
!           write(*,*) 'coucou',line
281
           do i=1,nnn
282
              read(ifil,'(A)',err=99,end=99) Line
283
              do while (line(1:1)==' ')
284
                 line=line(2:)
285
              end do
286
              atoms(i)=line(1:3)
287
!            write(*,*) 'coucou atoms',atoms(i)
288
              do while (line(1:1).ne.' ')
289
                 line=line(2:)
290
              end do
291
!              WRITE(*,*) 'coucou2:',i,line
292
              read(line,*) (tab(j,i),j=1,3)
293
           end do
294
!           WRITE(*,*) 'coucou:',nnn,tab(1,1)
295
           return
296
99           nnn=0
297
!             WRITE(*,*) 'Erreur lecture',ifil,nnn
298
             return
299
          end
300

    
301
!--------------------------------------------------------------
302

    
303
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
304
! READLINE
305
! This subroutine reads a line for a file, and converts
306
! the first field into a character variable, and the rest into 4 integers
307
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
308

    
309
      SUBROUTINE READLINE(IOIN,Type,Line)
310
        
311
       IMPLICIT NONE
312
      CHARACTER(5) :: Type
313
      CHARACTER(120) :: LINE
314
      INTEGER(4) :: i,IOIN
315
      READ(IOIN,'(A120)',ERR=999,END=999) LINE
316
!      WRITE(*,*) Line
317
      DO WHILE (LINE(1:1).eq.' ')
318
         LINE=LINE(2:)
319
      END DO
320

    
321
      i=1
322
      DO WHILE (LINE(i:i).ne.' ')
323
         i=i+1
324
      END DO
325
      
326
      if (i.ge.6) THEN
327
         WRITE(*,*) 'Pb with READLINE:',LINE
328
         GOTO 999
329
      END IF
330
      Type=LINE(1:i-1)
331
      LINE=LINE(i:120) // " 0 0 0 0"
332

    
333
      RETURN
334
 999  Type="E"
335
      END
336

    
337

    
338

    
339
      SUBROUTINE Analyse(geos)
340

    
341
      IMPLICIT NONE
342
          INCLUDE "Xyz2Path.param"
343

    
344
      REAL(8) :: geos(3,maxnat)
345
      REAL(8) :: AU2PS,Pi
346
      INTEGER(4) :: i,j,k
347
      REAL(8) :: V1x,V1y,V1z,V2x,V2y,V2z,V3x,V3y,V3z
348
      REAL(8) :: d1,d2,ca,sa
349
      REAL(8) :: V4x,v4y,v4z,v5x,v5y,v5z
350
      REAL(8) :: COG(3)
351
      LOGICAL :: Debug=.FALSE.
352
      INTEGER(4) :: Nat,NbDist, NbAngle, NbDie,NbCOM
353
      INTEGER(4) :: At1B(MaxNat),At2B(MaxNat)
354
      INTEGER(4) :: At1A(MaxNat),At2A(MaxNat),At3A(MaxNat)
355
      INTEGER(4) :: At1D(MaxNat),At2D(MaxNat),At3D(MaxNat),At4D(MaxNat)
356
      INTEGER(4) :: AtCom(0:MaxNat,MaxNat)
357
      REAL(8) ::    VB(MaxNat),VA(MaxNat),VD(MaxNat),VCOM(MaxNat)
358

    
359
      COMMON /Indices/Nat,NbDist,NbAngle,NbDie,NbCom,      &
360
           At1B,At2B,At1A,At2A,At3A,At1D,At2D,At3D,At4D,AtCom
361
      COMMON /Values/VB,VA,VD,VCom
362
      COMMON /Const/AU2ps,Pi
363

    
364
      if (debug) THEN
365
         DO I=1,Nat
366
            WRITE(*,'(1X,I3,3(1X,F15.6))') i,(geos(j,i),j=1,3)
367
         END DO
368
      END IF
369

    
370
! First, we create the Centre of Mass atoms
371
            DO i=1,NbCOM
372
              COG(1)=0.
373
              COG(2)=0.
374
              COG(3)=0.
375
              DO j=1,AtCOm(0,i)
376
                 DO k=1,3
377
                    COG(k)=COG(k)+geos(k,AtCom(j,i))
378
                 END DO
379
              END DO
380
              DO k=1,3
381
                 COG(k)=COG(k)/AtCOM(0,i)
382
                 geos(k,Nat+i)=COG(k)
383
              END DO
384
            END DO
385

    
386

    
387
            DO i=1,NbDist
388
               VB(i)=sqrt((geos(1,At1B(i))-geos(1,At2B(i)))**2+  &
389
           (geos(2,At1B(i))-geos(2,At2B(i)))**2+                &      
390
           (geos(3,At1B(i))-geos(3,At2B(i)))**2)
391
            END DO
392
            DO i=1,NbAngle
393
               v1x=geos(1,At1A(i))-geos(1,At2A(i))
394
               v1y=geos(2,At1A(i))-geos(2,At2A(i))
395
               v1z=geos(3,At1A(i))-geos(3,At2A(i))
396
               d1=sqrt(v1x**2+v1y**2+v1z**2)
397
               v2x=geos(1,At3A(i))-geos(1,At2A(i))
398
               v2y=geos(2,At3A(i))-geos(2,At2A(i))
399
               v2z=geos(3,At3A(i))-geos(3,At2A(i))
400
               d2=sqrt(v2x**2+v2y**2+v2z**2)
401
               VA(i)=acos((v1x*v2x+v1y*v2y+v1z*v2z)/(d1*d2))*180./Pi
402
            END DO
403
            DO i=1,NbDie
404
!               WRITE(*,*) At1D(i),At2D(i),At3D(i),At4D(i)
405
!              WRITE(*,*) geos(1,At1D(i)),geos(2,At1D(i)),geos(3,At1D(i))
406
!              WRITE(*,*) geos(1,At2D(i)),geos(2,At2D(i)),geos(3,At2D(i))
407
               v1x=geos(1,At1D(i))-geos(1,At2D(i))
408
               v1y=geos(2,At1D(i))-geos(2,At2D(i))
409
               v1z=geos(3,At1D(i))-geos(3,At2D(i))
410
               v2x=geos(1,At3D(i))-geos(1,At2D(i))
411
               v2y=geos(2,At3D(i))-geos(2,At2D(i))
412
               v2z=geos(3,At3D(i))-geos(3,At2D(i))
413
               v3x=geos(1,At4D(i))-geos(1,At3D(i))
414
               v3y=geos(2,At4D(i))-geos(2,At3D(i))
415
               v3z=geos(3,At4D(i))-geos(3,At3D(i))
416
          
417
               v4x=v1y*v2z-v1z*v2y
418
               v4y=v1z*v2x-v1x*v2z
419
               v4z=v1x*v2y-v1y*v2x
420
               d1=sqrt(v4x**2+v4y**2+v4z**2)
421
               v5x=-v2y*v3z+v2z*v3y
422
               v5y=-v2z*v3x+v2x*v3z
423
               v5z=-v2x*v3y+v2y*v3x
424
               d2=sqrt(v5x**2+v5y**2+v5z**2)
425
               ca=(v4x*v5x+v4y*v5y+v4z*v5z)/(d1*d2)
426
               sa=v1x*v5x+v1y*v5y+v1z*v5z
427
               VD(i)=acos(ca)*180./Pi
428
               if (sa<0.) VD(i)=-VD(i)
429
!               WRITE(*,*) "Dihe",v5x,v5y,v5z,v4x,v4y,v4z,d1,d2,
430
!     &(v4x*v5x+v4y*v5y+v4z*v5z)/(d1*d2),pi
431
!!!!!!!!! Another solution, more elegant ?
432
!               norm2=sqrt(v2x**2+v2y**2+v2z**2)
433
!               sa=(v4x*(v5y*v2z-v5z*v2y) 
434
!     *            -v4y*(v5x*v2z-v5z*v2x) 
435
!     *            +v4z*(v5x*v2y-v5y*v2x))
436
!     *       /(d1*norm2*d2)
437
!               angle_d=datan2(sa,ca)*180./Pi
438
!               WRITE(*,*) sa,ca,angle_d,d1,d2,norm2
439
!               WRITE(*,*) VD(i),angle_d
440
            END DO
441

    
442
          END SUBROUTINE Analyse
443

    
444
!C================================================================
445
!C    Convertit un nom d'atome (2 lettres) en nombre de masse (entier)
446
!C cette fonction a ete modifiee pour pouvoir convertir un nom
447
!C d'atome avec un numero a la fin...
448
!C================================================================
449

    
450
        FUNCTION ConvertNumAt(ATOM)
451

    
452

    
453
        IMPLICIT NONE 
454

    
455
        INTEGER(4) :: I,Long,ConvertNumAt,IC
456
        character(2) :: ATOM,L_Atom
457
        character(3) :: ATOME
458
        INTEGER(4), PARAMETER ::  Max_Z=86
459

    
460

    
461
  CHARACTER(2) :: Nom(0:max_Z)=(/ ' X',' H',                   'HE',  &
462
                'LI','BE',            ' B',' C',' N',' O',' F','NE',  &
463
                'NA','MG',            'AL','SI',' P',' S','CL','AR',  &
464
                ' K','CA',                                            &
465
            'SC','TI',' V','CR','MN','FE','CO','NI','CU','ZN',        &
466
                                      'GA','GE','AS','SE','BR','KR',  &
467
                'RB','SR',                                            &
468
            ' Y','ZR','NB','MO','TC','RU','RH','PD','AG','CD',        &
469
                                      'IN','SN','SB','TE',' I','XE',  &
470
                'CS','BA',                                            &
471
            'LA',                                                     &
472
              'CE','PR','ND','PM','SM','EU','GD','TB','DY','HO',      &
473
              'ER','TM','YB','LU',                                    &
474
                     'HF','TA',' W','RE','OS','IR','PT','AU','HG',    &
475
                                      'TL','PB','BI','PO','AT','RN'/)
476

    
477

    
478
!C Verifie qu'il n'y a que des lettres et des espaces dans ATOM
479
!        WRITE(*,*) 'DBG CNVNUMAT, ATOM:',ATOM
480

    
481
        L_atom=Atom
482
        IF (ATOM(1:1).LT.'A') L_ATOM(1:1)=' '
483
        IC=Ichar(ATOM(1:1))
484
        IF ((ic.le.123).AND.(ic.ge.97)) L_ATOM(1:1)=CHAr(IC-32)
485
        IF (ATOM(2:2).LT.'A') L_ATOM(2:2)=' '
486
        IC=Ichar(ATOM(2:2))
487
        IF ((ic.le.123).AND.(ic.ge.97)) L_ATOM(2:2)=CHAr(IC-32)
488
!C Justifie le nom sur la droite (et non sur la gauche comme souvent...)
489
         Long=INDEX(L_ATOM,' ')-1
490
         ATOME=' ' // L_ATOM
491
        IF (Long.EQ.1) L_ATOM=ATOME(1:2)
492
!        WRITE(*,*) 'DBG CNVNUMAT, L_ATOM:',L_ATOM
493
        I=max_Z
494
        DO WHILE ((nom(I).NE.L_ATOM) .AND. (I.GT.0))
495
         I=I-1
496
        END DO
497
        ConvertNumAT=I
498
        END
499
! This subroutine calculates RMSD using quaternions.
500
! It is based on the F90 routine bu E. Coutsias
501
! http://www.math.unm.edu/~vageli/homepage.html
502
! I (PFL) have just translated it, and I have changed the diagonalization
503
! subroutine.
504
! I also made some changes to make it suitable for Cart package.
505
!----------------------------------------------------------------------
506
!----------------------------------------------------------------------
507
! Copyright (C) 2004, 2005 Chaok Seok, Evangelos Coutsias and Ken Dill
508
!      UCSF, Univeristy of New Mexico, Seoul National University
509
! Witten by Chaok Seok and Evangelos Coutsias 2004.
510
                                                                                
511
! This library is free software; you can redistribute it and/or
512
! modify it under the terms of the GNU Lesser General Public
513
! License as published by the Free Software Foundation; either
514
! version 2.1 of the License, or (at your option) any later version.
515
!
516
                                                                                
517
! This library is distributed in the hope that it will be useful,
518
! but WITHOUT ANY WARRANTY; without even the implied warranty of
519
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
520
! Lesser General Public License for more details.
521
!
522
                                                                                
523
! You should have received a copy of the GNU Lesser General Public
524
! License along with this library; if not, write to the Free Software
525
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
526
!----------------------------------------------------------------------------
527

    
528
      subroutine CalcRmsd(na, geom,geom2,U,rmsd,FRot,FAlign,debug)
529
!-----------------------------------------------------------------------
530
!  This subroutine calculates the least square rmsd of two coordinate
531
!  sets coord1(3,n) and coord2(3,n) using a method based on quaternion.
532
!  It then calculate the  rotation matrix U and the centers of coord, and uses
533
! them to align the two molecules.
534
!-----------------------------------------------------------------------
535

    
536

    
537
      IMPLICIT NONE
538

    
539
      INCLUDE "Xyz2Path.param"
540

    
541
      INTEGER(4) :: IOIN, IOOUT, IOSCRN, IOPRNT, IOTAMP
542
      COMMON/IODEFS/IOIN,IOOUT,IOSCRN,IOPRNT,IOTAMP
543

    
544
      INTEGER(4) :: na
545
      REAL(8) ::  geom(3,MaxNAt), Geom2(3,MaxNAt)
546
      REAL(8) ::  U(3,3), rmsd
547
      LOGICAL  FRot,FAlign,Debug
548

    
549
      REAL(8) :: Coord1(3,MaxNAt), Coord2(3,MaxNAt)
550
      REAL(8) ::  xc1,yc1,zc1, xc2,yc2,zc2
551
      
552

    
553
      INTEGER(4) :: i, j,  ia
554
      REAL(8) :: x_norm, y_norm, lambda
555
      REAL(8) :: Rmatrix(3,3)
556
      REAL(8) :: S(4,4)
557
      REAL(8) :: EigVec(4,4), EigVal(4)
558

    
559

    
560

    
561
  ! calculate the barycenters, centroidal coordinates, and the norms
562
      x_norm = 0.0d0
563
      y_norm = 0.0d0
564
      xc1=0.
565
      yc1=0.
566
      zc1=0.
567
      xc2=0.
568
      yc2=0.
569
      zc2=0.
570
      do ia=1,na
571
         xc1=xc1+geom(1,ia)
572
         xc2=xc2+geom2(1,ia)
573
         yc1=yc1+geom(2,ia)
574
         yc2=yc2+geom2(2,ia)
575
         zc1=zc1+geom(3,ia)
576
         zc2=zc2+geom2(3,ia)
577
!         if (debug) WRITE(*,'(A,I5,4(1X,F10.4))') 'ia...',ia,x(ia),
578
!     &                        x2(ia),xc1,xc2
579
      END DO
580
      xc1=xc1/dble(na)
581
      yc1=yc1/dble(na)
582
      zc1=zc1/dble(na)
583
      xc2=xc2/dble(na)
584
      yc2=yc2/dble(na)
585
      zc2=zc2/dble(na)
586

    
587
      IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center1',xc1,yc1,zc1
588
      IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center2',xc2,yc2,zc2
589
      do i=1,na
590
         Coord1(1,i)=geom(1,i)-xc1
591
         Coord1(2,i)=geom(2,i)-yc1
592
         Coord1(3,i)=geom(3,i)-zc1
593
         Coord2(1,i)=geom2(1,i)-xc2
594
         Coord2(2,i)=geom2(2,i)-yc2
595
         Coord2(3,i)=geom2(3,i)-zc2
596
         x_norm=x_norm+Coord1(1,i)**2+Coord1(2,i)**2+Coord1(3,i)**2
597
         y_norm=y_norm+Coord2(1,i)**2+Coord2(2,i)**2+Coord2(3,i)**2
598
      end do
599

    
600
      IF (debug) THEN
601
         WRITE(*,*) "R matrix"
602
         DO I=1,3
603
            WRITE(*,*) (RMatrix(I,j),j=1,3)
604
         END DO
605
      END IF
606

    
607
  ! calculate the R matrix
608
      do i = 1, 3
609
         do j = 1, 3
610
            Rmatrix(i,j)=0.
611
            do ia=1,na
612
               Rmatrix(i,j) = Rmatrix(i,j)+Coord1(i,ia)*Coord2(j,ia)
613
            END DO
614
         end do
615
      end do
616

    
617
      IF (debug) THEN
618
         WRITE(*,*) "R matrix"
619
         DO I=1,3
620
            WRITE(*,*) (RMatrix(I,j),j=1,3)
621
         END DO
622
      END IF
623

    
624

    
625
  ! S matrix
626
      S(1, 1) = Rmatrix(1, 1) + Rmatrix(2, 2) + Rmatrix(3, 3)
627
      S(2, 1) = Rmatrix(2, 3) - Rmatrix(3, 2)
628
      S(3, 1) = Rmatrix(3, 1) - Rmatrix(1, 3)
629
      S(4, 1) = Rmatrix(1, 2) - Rmatrix(2, 1)
630

    
631
      S(1, 2) = S(2, 1)
632
      S(2, 2) = Rmatrix(1, 1) - Rmatrix(2, 2) - Rmatrix(3, 3)
633
      S(3, 2) = Rmatrix(1, 2) + Rmatrix(2, 1)
634
      S(4, 2) = Rmatrix(1, 3) + Rmatrix(3, 1)
635

    
636
      S(1, 3) = S(3, 1)
637
      S(2, 3) = S(3, 2)
638
      S(3, 3) =-Rmatrix(1, 1) + Rmatrix(2, 2) - Rmatrix(3, 3)
639
      S(4, 3) = Rmatrix(2, 3) + Rmatrix(3, 2)
640

    
641
      S(1, 4) = S(4, 1)
642
      S(2, 4) = S(4, 2)
643
      S(3, 4) = S(4, 3)
644
      S(4, 4) =-Rmatrix(1, 1) - Rmatrix(2, 2) + Rmatrix(3, 3) 
645

    
646

    
647
! PFL : I use my usual Jacobi diagonalisation
648
  ! Calculate eigenvalues and eigenvectors, and 
649
  ! take the maximum eigenvalue lambda and the corresponding eigenvector q.
650

    
651
      IF (debug) THEN
652
         WRITE(*,*) "S matrix"
653
         DO I=1,4
654
            WRITE(*,*) (S(I,j),j=1,4)
655
         END DO
656
      END IF
657

    
658
      Call Jacobi(S,4,EigVal,EigVec,4)
659

    
660
      Call Trie(4,EigVal,EigVec,4)
661

    
662
      Lambda=EigVal(4)
663

    
664
! RMS Deviation
665
       rmsd=sqrt(max(0.0d0,((x_norm+y_norm)-2.0d0*lambda))/dble(na))
666

    
667
      if (FRot.OR.FAlign) Call rotation_matrix(EigVec(1,4),U)
668
      IF (FAlign) THEN
669
         DO I=1,na
670
            geom2(1,i)=Coord2(1,i)*U(1,1)+Coord2(2,i)*U(2,1)+Coord2(3,i)*U(3,1) +xc1
671
            geom2(2,i)=Coord2(1,i)*U(1,2)+Coord2(2,i)*U(2,2)+Coord2(3,i)*U(3,2) +yc1
672
            geom2(3,i)=Coord2(1,i)*U(1,3)+Coord2(2,i)*U(2,3)+Coord2(3,i)*U(3,3) +zc1
673
         END DO
674
      END IF
675

    
676
      END
677

    
678

    
679
!-----------------------------------------------------------------------
680
      subroutine rotation_matrix(q, U)
681
!-----------------------------------------------------------------------
682
! This subroutine constructs rotation matrix U from quaternion q.
683
!-----------------------------------------------------------------------
684
! This subroutine calculates RMSD using quaternions.
685
! It is based on the F90 routine bu E. Coutsias
686
! http://www.math.unm.edu/~vageli/homepage.html
687
! I (PFL) have just translated it, and I have changed the diagonalization
688
! subroutine.
689
! I also made some changes to make it suitable for Cart package.
690
!----------------------------------------------------------------------
691
!----------------------------------------------------------------------
692
! Copyright (C) 2004, 2005 Chaok Seok, Evangelos Coutsias and Ken Dill
693
!      UCSF, Univeristy of New Mexico, Seoul National University
694
! Witten by Chaok Seok and Evangelos Coutsias 2004.
695
                                                                                
696
! This library is free software; you can redistribute it and/or
697
! modify it under the terms of the GNU Lesser General Public
698
! License as published by the Free Software Foundation; either
699
! version 2.1 of the License, or (at your option) any later version.
700
!
701
                                                                                
702
! This library is distributed in the hope that it will be useful,
703
! but WITHOUT ANY WARRANTY; without even the implied warranty of
704
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
705
! Lesser General Public License for more details.
706
!
707
                                                                                
708
! You should have received a copy of the GNU Lesser General Public
709
! License along with this library; if not, write to the Free Software
710
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
711
!----------------------------------------------------------------------------
712

    
713
      REAL(8) :: q(4)
714
      REAL(8) :: U(3,3)
715
      REAL(8) :: q0,q1,q2,q3,b0,b1,b2,b3,q00,q01,q02,q03
716
      REAL(8) :: q11,q12,q13,q22,q23,q33
717

    
718
      q0 = q(1)
719
      q1 = q(2)
720
      q2 = q(3)
721
      q3 = q(4)
722

    
723
      b0 = 2.0d0*q0
724
      b1 = 2.0d0*q1
725
      b2 = 2.0d0*q2
726
      b3 = 2.0d0*q3
727

    
728
      q00 = b0*q0-1.0d0
729
      q01 = b0*q1
730
      q02 = b0*q2
731
      q03 = b0*q3
732

    
733
      q11 = b1*q1
734
      q12 = b1*q2
735
      q13 = b1*q3  
736
      
737
      q22 = b2*q2
738
      q23 = b2*q3
739
      
740
      q33 = b3*q3 
741
      
742
      U(1,1) = q00+q11
743
      U(1,2) = q12-q03
744
      U(1,3) = q13+q02
745
      
746
      U(2,1) = q12+q03
747
      U(2,2) = q00+q22
748
      U(2,3) = q23-q01
749
      
750
      U(3,1) = q13-q02
751
      U(3,2) = q23+q01
752
      U(3,3) = q00+q33
753
      
754
      end 
755
      
756
!c============================================================
757
!c
758
!c ++  diagonalisation par jacobi
759
!c     vecteur propre i : V(i,i)
760
!c     valeur propre i : D(i)
761
!c
762
!c============================================================
763
!c 
764
      SUBROUTINE JACOBI(A,N,D,V,max_N)
765
      IMPLICIT REAL(8) (A-H,O-Z)
766
      parameter (max_it=500)
767
      DIMENSION A(max_N,max_N),B(max_N),Z(max_N)
768
      DIMENSION V(max_N,max_N),D(max_N)
769

    
770

    
771
      DO 12 IP=1,N
772
        DO 11 IQ=1,N
773
          V(IP,IQ)=0.
774
11      CONTINUE
775
        V(IP,IP)=1.
776
12    CONTINUE
777
      DO 13 IP=1,N
778
        B(IP)=A(IP,IP)
779
        D(IP)=B(IP)
780
        Z(IP)=0.
781
13    CONTINUE
782
      NROT=0
783
      DO 24 I=1,max_it
784
        SM=0.
785
        DO 15 IP=1,N-1
786
          DO 14 IQ=IP+1,N
787
            SM=SM+ABS(A(IP,IQ))
788
14        CONTINUE
789
15      CONTINUE
790
        IF(SM.EQ.0.)RETURN
791
        IF(I.LT.4)THEN
792
          TRESH=0.2*SM/N**2
793
        ELSE
794
          TRESH=0.
795
        ENDIF
796
        DO 22 IP=1,N-1
797
          DO 21 IQ=IP+1,N
798
            G=100.*ABS(A(IP,IQ))
799
            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) &
800
              .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
801
              A(IP,IQ)=0.
802
            ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
803
              H=D(IQ)-D(IP)
804
              IF(ABS(H)+G.EQ.ABS(H))THEN
805
                T=A(IP,IQ)/H
806
              ELSE
807
                THETA=0.5*H/A(IP,IQ)
808
                T=1./(ABS(THETA)+SQRT(1.+THETA**2))
809
                IF(THETA.LT.0.)T=-T
810
              ENDIF
811
              C=1./SQRT(1+T**2)
812
              S=T*C
813
              TAU=S/(1.+C)
814
              H=T*A(IP,IQ)
815
              Z(IP)=Z(IP)-H
816
              Z(IQ)=Z(IQ)+H
817
              D(IP)=D(IP)-H
818
              D(IQ)=D(IQ)+H
819
              A(IP,IQ)=0.
820
              DO 16 J=1,IP-1
821
                G=A(J,IP)
822
                H=A(J,IQ)
823
                A(J,IP)=G-S*(H+G*TAU)
824
                A(J,IQ)=H+S*(G-H*TAU)
825
16            CONTINUE
826
              DO 17 J=IP+1,IQ-1
827
                G=A(IP,J)
828
                H=A(J,IQ)
829
                A(IP,J)=G-S*(H+G*TAU)
830
                A(J,IQ)=H+S*(G-H*TAU)
831
17            CONTINUE
832
              DO 18 J=IQ+1,N
833
                G=A(IP,J)
834
                H=A(IQ,J)
835
                A(IP,J)=G-S*(H+G*TAU)
836
                A(IQ,J)=H+S*(G-H*TAU)
837
18            CONTINUE
838
              DO 19 J=1,N
839
                G=V(J,IP)
840
                H=V(J,IQ)
841
                V(J,IP)=G-S*(H+G*TAU)
842
                V(J,IQ)=H+S*(G-H*TAU)
843
19            CONTINUE
844
              NROT=NROT+1
845
            ENDIF
846
21        CONTINUE
847
22      CONTINUE
848
        DO 23 IP=1,N
849
          B(IP)=B(IP)+Z(IP)
850
          D(IP)=B(IP)
851
          Z(IP)=0.
852
23      CONTINUE
853
24    CONTINUE
854
      write(6,*) max_it,' iterations should never happen'
855
      STOP
856
      RETURN
857
      END
858
!
859
!============================================================
860
!c
861
!c ++  trie des vecteur dans l'ordre croissant
862
!c
863
!c============================================================
864
!c
865
        SUBROUTINE trie(nb_niv,ene,psi,max_niv)
866
        integer(4) :: i,j,k,nb_niv,max_niv
867
        REAL(8) :: ene(max_niv),psi(max_niv,max_niv)
868
        REAL(8) :: a
869

    
870
        DO i=1,nb_niv
871
          DO j=i+1,nb_niv
872
            IF (ene(i) .GT. ene(j)) THEN
873
!c              permutation
874
              a=ene(i)
875
              ene(i)=ene(j)
876
              ene(j)=a
877
              DO k=1,nb_niv
878
                 a=psi(k,i)
879
                 psi(k,i)=psi(k,j)
880
                 psi(k,j)=a
881
              END DO
882
            END IF
883
          END DO
884
        END DO
885

    
886
      END SUBROUTINE trie
887