Revision 5 src/Path.f90

Path.f90 (revision 5)
103 103
! When using internal coordinates (zmat, hybrid or mixes), Path calculates the number
104 104
! of fragments for each geometry. The reference one (ie the one used to compute the internal
105 105
! coordinates) is now the one with the least fragments.
106
! user can also choose one geometry by specifying IGeomRef in the Path namelis.
106
! user can also choose one geometry by specifying IGeomRef in the Path namelist.
107 107
! v3.96
108 108
! I have added a test for the geometry step: it does not allow valence angle to 
109 109
! be negative or bigger than 180. I hope that this might help geometry optimization
......
204 204
! decomposition, and to update this decomposition rather than the Hessian.
205 205
! See Nocedal&Wright p141.
206 206
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
207
!
207 208
! Project Name has changed to Opt'n Path
208 209
! OpenPaht 2011 May-June -> Opt'n Path v1.4 
210
!
211
! 2013 Jan 18 Opt'n Path v1.47 (going to 1.5 slowly)
212
! We add Siesta (and start thinking of add CP2K)
213
!  as a new "energy engine" code
214
! TO DO to go to 1.5: finish cleaning...
215
! 2013 Feb Opt'n Path v1.48
216
! We add some keyword for more output for Geometry Optimizations
217
! GeomFile: name of the file to print out the geometries and their energy
218
! as Xyz file. (only format for now)
209 219

  
210

  
211 220
PROGRAM PathOpt
212 221

  
213 222
  use VarTypes
......
216 225

  
217 226
  IMPLICIT NONE
218 227

  
219
  CHARACTER(132) :: FileIn, FileOut, Line,Line2
220
  CHARACTER(32) :: Input,poscar
228
  CHARACTER(132) :: FileIn, FileOut, Line
221 229

  
222 230
  INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF
223 231
  INTEGER(KINT) :: IOpt
224
  INTEGER(KINT) :: Idx, LineL,LenLine
232
  INTEGER(KINT) :: Idx
225 233
  INTEGER(KINT) :: N3at, NFree, Nfr
226
  ! Temporary integer counter
227
  INTEGER(KINT) :: NTmp
228 234

  
229 235
  INTEGER(KINT) :: List(2*MaxFroz)
230 236

  
231
  REAL(KREAL) :: E,FactStep,B(3) !E=Energy
237
  REAL(KREAL) :: E,FactStep !E=Energy
232 238
  REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad
233 239
  REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E
234 240
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord
......
248 254
  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF, where it is deallocated: Prakash 
249 255

  
250 256
! these are used to read temporary coordinates
251
  REAL(KREAL) :: xtmp,ytmp,ztmp
252

  
253 257
  LOGICAL :: FFrozen,FCart
254 258

  
255 259
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
......
262 266
  LOGICAL, ALLOCATABLE :: ListAtFroz(:)
263 267
  INTEGER(KINT) :: Iat1, Iat2, Iat3
264 268

  
265
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
266
  !
267
  ! For VASP
268
  !
269
  INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na
270
  INTEGER(KINT) :: NbType, NbTypeUser
271

  
272 269
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
273 270

  
274
  LOGICAL :: TChk
275 271
  LOGICAL :: Debug, Fini,Flag_Opt_Geom
276 272

  
277 273
!  INTEGER(KINT), EXTERNAL :: Iargc
......
295 291
       CHARACTER(32), INTENT(IN) :: input
296 292
     END SUBROUTINE READ_GEOM
297 293

  
294

  
298 295
     subroutine hupdate_all (n, dx, dgrad, HessOld)
299 296

  
300 297
       IMPLICIT NONE
......
330 327
       OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, &
331 328
       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &
332 329
       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &
333
       NGintMax, Align, CalcEReac,CalcEProd
330
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile
334 331

  
335 332
  Namelist/cartlist/list
336 333
  Namelist/frozenlist/list
......
342 339
  ! 
343 340
  ! We print the Version info in file
344 341
  !
345
  WRITE(*,'(1X,A)') "OpenPath v1.4 (c) PFL/PD 2007-2011"
342
  WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 2007-2013"
346 343

  
347 344
  ! We read the files names
348 345
!  SELECT CASE (iargc())
......
371 368
        WRITE(*,*) "Input: string that indicates the type of the input geometries."
372 369
        WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp"
373 370
        WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
374
        WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac or Ext"
371
        WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac, Test, Siesta or Ext"
375 372
        WRITE(*,*) "      In case of a Gaussian or Mopac calculations, input must be set to Cart."
376 373
        WRITE(*,*) "      One example of a gaussian/mopac input should be added at the end of the"
377 374
        WRITE(*,*) "      input file.See example file Test_G03.path or Test_Mopac.path"
......
409 406
!        WRITE(*,*) " HesUpd: Deprecated. Use Hupdate."
410 407
        WRITE(*,*) " HUpdate: method to update the hessian. By default, it is Murtagh-Sargent"
411 408
        WRITE(*,*) "      Except for geometry optimization where it is BFGS."
409
        WRITE(*,*) " GeomFile: Name of the file to print the geometries and their energy"
410
        WRITE(*,*) "      during a geometry optimization. If this variable is not given"
411
        WRITE(*,*) "      then nothing is printed."
412 412
        WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"
413 413
        WRITE(*,*) "Smax: Maximum length of a step during path optimization"
414 414
        WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
......
461 461
        WRITE(*,*) ""
462 462
        STOP
463 463
     END IF
464

  
464 465
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
465 466
     IOOUT=6
466 467
  CASE (0)
......
496 497
  IniHUp=.FALSE.
497 498
  HupNeighbour=.FALSE.
498 499
  HesUpd="EMPTY"
499
  HUpdate="MS"
500
  HUpdate="EMPTY"
500 501
  FFrozen=.FALSE.
501 502
  FCart=.FALSE.
502 503
  List=0
......
508 509
  EReac=0.d0
509 510
  EProd=0.d0
510 511
  OptGeom=-1
512
  GeomFile="EMPTY"
511 513
  PathOnly=.FALSE.
512 514
  Prog='EMPTY'
513 515
  ProgExe='EMPTY'
......
532 534
  ! Read the namelist
533 535
  read (IOIN,path)
534 536

  
535
  debug=valid("Main")
537
  debug=valid("Main").or.valid("Path")
538

  
536 539
  FCalcHess=valid("CalcHess")
537 540
  PathName=AdjustL(PathName)
538 541
  Coord=AdjustL(Coord)
......
569 572
     STOP
570 573
  END IF
571 574

  
575
  IF (HUpdate=='EMPTY') THEN
576
     IF (OptGeom>=1) THEN
577
        HupDate="BFGS"
578
     ELSE
579
        HUpdate="MS"
580
     END IF
581
  END IF
572 582

  
583
  If (COORD.EQ.'XYZ') THEN 
584
     COORD='CART'
585
  END IF
586

  
573 587
  IF (COORD.EQ.'CART') THEN
574 588
     Renum=.FALSE.
575
     IGeomRef=1
589
     IF (OptGeom>0) THEN
590
        IGeomRef=OptGeom
591
     ELSE
592
        IGeomRef=1
593
     END IF
576 594
  END IF
577 595

  
578
  if (Prog.EQ.'EMPTY') THEN
579
     Prog='GAUSSIAN'
580
     If (ProgExe=="EMPTY") ProgExe="g09"
581
  END IF
582
  if (Prog.EQ.'G03') THEN
583
     Prog='GAUSSIAN'
584
     If (ProgExe=="EMPTY") ProgExe="g03"
585
  END IF
586
  if (Prog.EQ.'G09') THEN
587
     Prog='GAUSSIAN'
588
     If (ProgExe=="EMPTY") ProgExe="g09"
589
  END IF
590

  
591
  if (Prog.EQ.'TM') Prog="TURBOMOLE"
596
  
592 597
  if (Prog.EQ.'TEST2') Prog="CHAMFRE"
598
  if (Prog.EQ.'2D') Prog="TEST2D"
593 599
  if (Prog.EQ.'TEST3') Prog="LEPS"
594 600

  
595 601

  
596 602
  Select Case (Prog)
603
    CASE ("EMPTY","G09","GAUSSIAN")
604
     Prog='GAUSSIAN'
605
     If (ProgExe=="EMPTY") ProgExe="g09"
606
     if (CalcName=="EMPTY") CalcName="Path"
607
     if (ISuffix=="EMPTY")  ISuffix="com"
608
     if (OSuffix=="EMPTY")  OSuffix="log"
609
    CASE ("G03")
610
     Prog='GAUSSIAN'
611
     If (ProgExe=="EMPTY") ProgExe="g03"
612
     if (CalcName=="EMPTY") CalcName="Path"
613
     if (ISuffix=="EMPTY")  ISuffix="com"
614
     if (OSuffix=="EMPTY")  OSuffix="log"
597 615
    CASE ("EXT")
616
     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"
598 617
     if (CalcName=="EMPTY") CalcName="Ext"
599 618
     if (ISuffix=="EMPTY")  ISuffix="in"
600 619
     if (OSuffix=="EMPTY")  OSuffix="out"
601 620
    CASE ('MOPAC')
621
     If (ProgExe=="EMPTY") ProgExe="mopac"
602 622
     if (CalcName=="EMPTY") CalcName="Path"
603 623
     if (ISuffix=="EMPTY")  ISuffix="mop"
604 624
     if (OSuffix=="EMPTY")  OSuffix="out"
605
    CASE ("GAUSSIAN")
625
    CASE ("SIESTA")
626
     If (ProgExe=="EMPTY") ProgExe="siesta"
606 627
     if (CalcName=="EMPTY") CalcName="Path"
607
     if (ISuffix=="EMPTY")  ISuffix="com"
608
     if (OSuffix=="EMPTY")  OSuffix="log"
609
    CASE DEFAULT
610
     if (CalcName=="EMPTY") CalcName="Path"
611
     if (ISuffix=="EMPTY")  ISuffix="in"
628
     if (ISuffix=="EMPTY")  ISuffix="fdf"
612 629
     if (OSuffix=="EMPTY")  OSuffix="out"
630
    CASE ("VASP")
631
     If (ProgExe=="EMPTY") ProgExe="VASP"
632
    CASE ("TM","TURBOMOLE")
633
     Prog="TURBOMOLE"
634
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
635
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
636
     ProgExe=""
637
    CASE DEFAULT
638
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
639
     STOP
613 640
  END SELECT 
614 641

  
615 642
  if (Input.EQ.'EMPTY') THEN
......
624 651
    WRITE(*,*) "Input has been set to the default: ",INPUT
625 652
 END IF
626 653

  
654

  
655

  
627 656
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
628 657

  
629 658
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
......
632 661
  END IF
633 662

  
634 663
  if (OptGeom.GE.1) THEN
635
     HUpdate="BFGS"
664
     FPrintGeom=.FALSE.
665
     If (GeomFile/='EMPTY') THEN
666
        FPrintGeom=.TRUE.
667
     END IF
636 668
     IF (IGeomRef.LE.0) THEN
637 669
        IGeomRef=OptGeom
638 670
     ELSE
......
647 679
     END IF
648 680
  END IF
649 681

  
650
  IF (ProgExe.EQ.'EMPTY') THEN
651
     SELECT CASE (PROG)
652
     CASE ("GAUSSIAN")
653
        ProgExe="g03"
654
     CASE ("MOPAC")
655
        ProgExe="mopac"
656
     CASE ("EXT")
657
        ProgExe="./Ext.exe"
658
     CASE ("VASP")
659
        ProgExe='VASP'
660
     CASE ("TURBOMOLE")
661
        ProgExe='jobex -c 1 -keep'
662
     CASE ("TEST")
663
        ProgExe=""
664
     CASE ("CHAMFRE")
665
        ProgExe=""
666
     CASE DEFAULT
667
        WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
668
        STOP
669
     END SELECT
670
  END IF
682
  ALLOCATE(Cart(Nat))
683
  Cart=0
671 684

  
672
  if (FCart.AND.FFrozen) THEN
673
     ! We have to be carreful because user might use any order
674
     ! to give the two lists CartList and FrozenList
675
     READ(IOIN,'(A)') Line
676
     Call Upcase(Line)
677
     if (Index(Line,"CART").NE.0) THEN
678
        ALLOCATE(Cart(Nat))
679
        BackSpace(IOIN)
680
        List=0
681
        READ(IOIN,cartlist)
682
        IF (COORD.NE.'CART') THEN
683
           Cart=List(1:Nat)
684
           if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
685
        ELSE
686
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
687
           WRITE(*,*) "I will discard the cartlist"
688
           Cart=0
689
        END IF
685
  ALLOCATE(Frozen(Nat))
686
  Frozen=0
690 687

  
691
        ALLOCATE(Frozen(Nat))
692
        BACKSPACE(IOIN)
693
        List=0
694
        READ(IOIN,Frozenlist)
695
        Frozen=List(1:Nat)
696
        if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
688
  IF (FCart) THEN
689
! We rewind the file to be sure to browse all of it to read the Cart namelist
690
     REWIND(IOIN)
691
     List=0
692
     READ(IOIN,cartlist)
693
     IF (COORD.NE.'CART') THEN
694
        Cart=List(1:Nat)
695
        if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
697 696
     ELSE
698
        ALLOCATE(Frozen(Nat))
699
        BACKSPACE(IOIN)
700
        List=0
701
        READ(IOIN,Frozenlist)
702
        Frozen=List(1:Nat)
703
        if (debug) WRITE(*,*) "DBG Path L523 Frozen",Frozen
704
        ALLOCATE(Cart(Nat))
705
        BackSpace(IOIN)
706
        List=0
707
        READ(IOIN,cartlist)
708
        IF (COORD.NE.'CART') THEN
709
           Cart=List(1:Nat)
710
           if (debug) WRITE(*,*) "DBG Path L548 Cart",Cart
711
        ELSE
712
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
713
           WRITE(*,*) "I will discard the cartlist"
714
           Cart=0
715
        END IF
716
     END IF
717
  ELSE
718
     if (FCart) THEN
719
        ALLOCATE(Cart(Nat))
720
        List=0
721
        READ(IOIN,cartlist)
722
        IF (COORD.NE.'CART') THEN
723
           Cart=List(1:Nat)
724
           if (debug) WRITE(*,*) "DBG Path L561 Cart",Cart
725
        ELSE
726
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
727
           WRITE(*,*) "I will discard the cartlist"
728
           Cart=0
729
        END IF
730
     ELSE
731
        IF ((PROG=="VASP").OR.(AutoCart)) THEN
732
           ALLOCATE(Cart(Nat))
733
        ELSE
734
           ALLOCATE(Cart(1))
735
        END IF
697
        WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
698
        WRITE(*,*) "I will discard the cartlist"
736 699
        Cart=0
737 700
     END IF
738

  
739
     if (FFrozen) THEN
740
        ALLOCATE(Frozen(Nat))
741
        List=0
742
        READ(IOIN,Frozenlist)
743
        Frozen=List(1:Nat)
744
        WRITE(*,*) Frozen
745
        if (debug) WRITE(*,*) "DBG Path L549 Frozen",Frozen
746
     ELSE
747
        IF (PROG=="VASP") THEN
748
           ALLOCATE(Frozen(Nat))
749
        ELSE
750
           ALLOCATE(Frozen(1))
751
        END IF
752
        Frozen=0
753
     END IF
754 701
  END IF
755 702

  
703
  if (FFrozen) THEN
756 704

  
705
     REWIND(IOIN)
706
     List=0
707
     READ(IOIN,Frozenlist)
708
     Frozen=List(1:Nat)
709
     if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
710
  END IF
711

  
757 712
  If (FCart) Call Expand(Cart,NCart,Nat)
758 713
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)
759 714

  
......
806 761
     Atome(I)=ConvertNumAt(AtName(I))
807 762
  END DO
808 763

  
764
!*********** HERE *************
765
!* To DO: 
766
!* deplacer le test sur les frozen dans une sous routine
767
!* deplacer la lecture de sinput dans une sous routine
768
!* creer lecture input siesta
769
!* creer egrad siesta
770
!*
771
!*************************************
809 772

  
810 773
  ! PFL 23 Sep 2008 ->
811 774
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
......
933 896
  N3at=3*Nat
934 897
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.
935 898

  
936
  ! if prog=gaussian, we have to read a full input
937
  ! to mimic it later
938
  if (prog=='GAUSSIAN') THEN
939
     ! We read the Gaussian input file
940
     ! First, the root
941
     IF (DEBUG) WRITE(*,*) "Reading Gauss Root"
942
     ALLOCATE(Gauss_Root)
943
     NULLIFY(Gauss_Root%next)
944
     Current => Gauss_root
945
     LineL=1
946
     DO WHILE (LineL.NE.0)
947
        READ(IOIN,'(A)') Line
948
        Line=AdjustL(Line)
949
        LineL=len(Trim(Line))
950
        Idx=INDEX(Line,"chk")
951
        IF ((LineL.NE.0).AND.(Idx.EQ.0)) THEN
952
           current%Line=TRIM(Line)
953
           ALLOCATE(current%next)
954
           Current => Current%next
955
           Nullify(Current%next)
956
        END IF
957
     END DO
958 899

  
959
!     Current => Gauss_root
960
!     DO WHILE (ASSOCIATED(Current%next))
961
!        WRITE(*,'(1X,A)') Trim(current%line)
962
!        Current => current%next
963
!     END DO
900
  Call ReadInput
964 901

  
965
     ! Now the comment... 
966
     IF (DEBUG) WRITE(*,*) "Reading Gauss Comment"
967
     ALLOCATE(Gauss_Comment)
968
     NuLLIFY(Gauss_Comment%Next)
969
     Current => Gauss_comment
970
     LineL=1
971
     DO WHILE (LineL.NE.0)
972
        READ(IOIN,'(A)') Line
973
        Line=AdjustL(Line)
974
        LineL=len(Trim(Line))
975
        IF (LineL.NE.0) THEN
976
           current%Line=TRIM(Line)
977
           ALLOCATE(current%next)
978
           Current => Current%next
979
           Nullify(Current%next)
980
        END IF
981
     END DO
982

  
983
 !    Current => Gauss_comment
984
 !    DO WHILE (ASSOCIATED(Current%next))
985
 !       WRITE(*,'(1X,A)') Trim(current%line)
986
 !       Current => current%next
987
 !    END DO
988

  
989
     ! Now the charge
990
     IF (DEBUG) WRITE(*,*) "Reading Gauss Charge"
991
     READ(IOIN,'(A)') Gauss_Charge
992
!     WRITE(*,*) Gauss_charge
993
     ! We now read the Paste part...
994
     ALLOCATE(Gauss_Paste(NAt))
995
     LenLine=1
996
     Iat=0
997
     Gauss_paste=" "
998
     DO While (LenLine.GT.0)
999
        READ(IOIN,'(A)') Line
1000
        Line=AdjustL(Line)
1001
        LenLine=Len_TRIM(Line)
1002
        IF (LenLine.GT.0) Iat=Iat+1
1003
        Idx=Index(Line,'.',BACK=.TRUE.)
1004
        !          WRITE(*,*) Idx,TRIM(Line)
1005
        !          WRITE(*,*) Idx,TRIM(Line(Idx+1:))
1006
        Line=ADJUSTL(Line(Idx+1:))
1007
        Idx=Index(Line," ")
1008
        LineL=LEN_TRIM(Line(Idx:))
1009
        !          WRITE(*,*) LineL,TRIM(Line(Idx:))
1010
        IF (LineL.GT.0) THEN 
1011
           Gauss_paste(Iat)=Line(Idx:)
1012
        END IF
1013
     END DO
1014
     IF (Iat.NE.Nat) THEN
1015
        WRITE(*,*) "I found ", Iat," lines for the geometry instead of ",Nat
1016
        WRITE(*,*) "ERROR. STOP"
1017
        WRITE(IOOUT,*) "I found ", Iat," lines for the geometry instead of ",Nat
1018
        WRITE(IOOUT,*) "ERROR. STOP"
1019
        STOP
1020
     END IF
1021
     ! We now read the last part
1022
     IF (DEBUG) WRITE(*,*) "Reading Gauss End"
1023
     !     READ(IOIN,'(A)') Line
1024
     ALLOCATE(Gauss_End)
1025
     NuLLIFY(Gauss_End%Next)
1026
     Current => Gauss_End
1027
     LineL=1
1028
     DO WHILE (1.EQ.1)
1029
        READ(IOIN,'(A)',END=999) Line
1030
        Line=AdjustL(Line)
1031
        LineL=len(Trim(Line))
1032
        current%Line=TRIM(Line)
1033
        ALLOCATE(current%next)
1034
        Current => Current%next
1035
        Nullify(Current%next)
1036
     END DO
1037
999  CONTINUE
1038

  
1039
     !     ! Write the gaussian input file for testing purposes
1040
     !     Current => Gauss_root
1041
     !     DO WHILE (ASSOCIATED(Current%next))
1042
     !        WRITE(*,'(1X,A)') Trim(current%line)
1043
     !        Current => current%next
1044
     !     END DO
1045

  
1046
     !        WRITE(*,*) 
1047
     !        WRITE(*,*) 'Comment original'
1048

  
1049
     !        Current => Gauss_comment
1050
     !        DO WHILE (ASSOCIATED(Current%next))
1051
     !           WRITE(*,'(1X,A)') Trim(current%line)
1052
     !           Current => current%next
1053
     !        END DO
1054

  
1055
     !        WRITE(*,*) 
1056
     !        WRITE(*,*) Trim(Gauss_charge)
1057

  
1058
     !        DO I=1,Nat
1059
     !           WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XyzGeomI(1:3,I,1), TRIM(Gauss_Paste(I))
1060
     !        END DO
1061

  
1062
     !        WRITE(*,*) 
1063
     !        Current => Gauss_End
1064
     !        DO WHILE (ASSOCIATED(Current%next))
1065
     !           WRITE(*,'(1X,A)') Trim(current%line)
1066
     !           Current => current%next
1067
     !        END DO
1068

  
1069
     !        WRITE(*,*) 
1070

  
1071
  END IF
1072

  
1073
! if prog=mopac, we have to read a full input to reproduce it later
1074
! The structure is:
1075
! A MOPAC data set normally consists of one line of keywords, two lines of user-defined text, then the coordinates
1076
! Then  a blank line or a line of 0.
1077
! then the symmetry description.
1078
! comment lines start with * and can be anywhere !!!
1079

  
1080
  if (prog=='MOPAC') THEN
1081
     ! We read the Gaussian input file
1082
     ! First, the root
1083
     IF (DEBUG) WRITE(*,*) "Reading Mopac input"
1084
     ALLOCATE(Mopac_Root)
1085
     NULLIFY(Mopac_Root%next)
1086
     ALLOCATE(Mopac_Comment)
1087
     NULLIFY(Mopac_Comment%next)
1088
     ALLOCATE(Mopac_End)
1089
     NuLLIFY(Mopac_End%Next)
1090
     Current => Mopac_root
1091
     CurCom => Mopac_Comment
1092
     LineL=1
1093
     NTmp=0
1094
     DO WHILE (NTmp.LT.3)
1095
        READ(IOIN,'(A)') Line
1096
        Line=AdjustL(Line)
1097
        LineL=len(Trim(Line))
1098
        IF (Line(1:1)/="*") THEN
1099
           IF (NTmp==0) THEN
1100
              Line2=Line
1101
              Call UpCase(Line2)
1102
              Idx=Index(Line2,'GRADIENTS')
1103
              If (Idx==0) Line=TRIM(Line) // " GRADIENTS"
1104
              Idx=Index(Line2,'1SCF')
1105
              If (Idx==0) Line=TRIM(Line) // " 1SCF"
1106
           END IF
1107
           current%Line=TRIM(Line)
1108
           ALLOCATE(current%next)
1109
           Current => Current%next
1110
           Nullify(Current%next)
1111
           NTmp=NTmp+1
1112
        ELSE
1113
           CurCom%Line=TRIM(LINE)
1114
           ALLOCATE(CurCom%Next)
1115
           CurCom => CurCom%Next
1116
           NULLIFY(CurCom%Next)
1117
        END IF
1118
     END DO
1119

  
1120
!     Current => Mopac_root
1121
!     DO WHILE (ASSOCIATED(Current%next))
1122
!        WRITE(*,'(1X,A)') Trim(current%line)
1123
!        Current => current%next
1124
!     END DO
1125

  
1126
! Now the geometry... that we just skip :)
1127
     IF (DEBUG) WRITE(*,*) "Reading Mopac Geometry"
1128
     Mopac_EndGeom=""
1129
     LineL=1
1130
     DO WHILE (LineL.NE.0)
1131
        READ(IOIN,'(A)',END=989) Line
1132
        Line=AdjustL(Line)
1133
        LineL=len(Trim(Line))
1134
! The last line might be either blank or filled with 0
1135
        IF (Line(1:1)=="0") THEN
1136
           LineL=0
1137
           Mopac_EndGeom=Trim(Line)
1138
        END IF
1139
        IF (Line(1:1)=="*") THEN
1140
           CurCom%Line=TRIM(LINE)
1141
           ALLOCATE(CurCom%Next)
1142
           CurCom => CurCom%Next
1143
           NULLIFY(CurCom%Next)
1144
        END IF
1145
     END DO
1146

  
1147
! If we are here, there might be something else to read: Mopac_end
1148

  
1149
     ! We now read the last part
1150
     IF (DEBUG) WRITE(*,*) "Reading Mopac End"
1151
     !     READ(IOIN,'(A)') Line
1152
     Current => Mopac_End
1153
     LineL=1
1154
     DO WHILE (1.EQ.1)
1155
        READ(IOIN,'(A)',END=989) Line
1156
        Line=AdjustL(Line)
1157
        LineL=len(Trim(Line))
1158
        IF (Line(1:1)/="*") THEN
1159
           current%Line=TRIM(Line)
1160
           ALLOCATE(current%next)
1161
           Current => Current%next
1162
           Nullify(Current%next)
1163
           NTmp=NTmp+1
1164
        ELSE
1165
           CurCom%Line=TRIM(LINE)
1166
           ALLOCATE(CurCom%Next)
1167
           CurCom => CurCom%Next
1168
           NULLIFY(CurCom%Next)
1169
        END IF
1170
     END DO
1171
989  CONTINUE
1172

  
1173
  END IF ! IF (PROG="MOPAC")
1174

  
1175
  if ((Prog=="VASP").AND.(input/="VASP")) THEN
1176
     ! Input was not Vasp, so many parameters are missing like lattice 
1177
     ! constants...
1178
     ! we read them now !
1179
     ALLOCATE(FFF(3,nat))
1180
     ! First geometry is a bit special for VASP as we have to set
1181
     ! many things
1182
     IF (DEBUG) WRITE(*,*) "Reading Vasp Parameters"
1183
     READ(IOIN,'(A)') Vasp_Title
1184
     READ(IOIN,*) Vasp_param
1185

  
1186
     READ(IOIN,*) Lat_a
1187
     READ(IOIN,*) Lat_b
1188
     READ(IOIN,*) Lat_c
1189

  
1190
     Lat_a=Lat_a*Vasp_param
1191
     Lat_b=Lat_b*Vasp_param
1192
     Lat_c=Lat_c*Vasp_param
1193

  
1194
     ALLOCATE(NbAtType(nat))
1195
     READ(IOIN,'(A)') Vasp_types
1196
     ! First, how many different types ?
1197
     NbAtType=0
1198
     READ(Vasp_types,*,END=998) NbAtType
1199
998  CONTINUE
1200
     NbType=0
1201
     DO WHILE (NbAtType(NbType+1).NE.0)
1202
        NbType=NbType+1
1203
     END DO
1204

  
1205
     ! Do we know the atom types ?
1206
     IF (AtTypes(1).EQ.'  ') THEN
1207
        ! user has not provided atom types... we have to find them ourselves
1208
        ! by looking into the POTCAR file...
1209
        INQUIRE(File="POTCAR", EXIST=Tchk)
1210
        IF (.NOT.Tchk) THEN
1211
           WRITE(*,*) "ERROR! No AtTypes provided, and POTCAR file not found"
1212
           STOP
1213
        END IF
1214
        OPEN(IOTMP,File="POTCAR")
1215
        DO I=1,NbType
1216
           Line='Empty'
1217
           DO WHILE (Line(1:2).NE.'US')
1218
              READ(IOTMP,'(A)') Line
1219
              Line=AdjustL(Line)
1220
           END DO
1221
           Line=adjustl(Line(3:))
1222
           AtTypes(I)=Line(1:2)
1223
        END DO
1224
        if (debug) WRITE(*,'(A,100(1X,A2))') "ReadG:VASP AtTypes",AtTypes(1:NbType)
1225
        CLOSE(IOTMP)
1226

  
1227
     ELSE  !AtTypes(1).EQ.'  '
1228
        ! user has provided atom types
1229
        NbTypeUser=0
1230
        DO WHILE (AtTypes(NbTypeUser+1).NE.'  ')
1231
           NbTypeUser=NbTypeUser+1
1232
        END DO
1233
        IF (NbType.NE.NbTypeUser) THEN
1234
           WRITE(*,*) "ERROR Read_Geom : NbType in POSCAR do not match AtTypes"
1235
           STOP
1236
        END IF
1237
     END IF
1238

  
1239
     IAt=1
1240
     DO I=1,NbType
1241
        DO J=1,NbAtType(I)
1242
           AtName(Iat)=AtTypes(I)
1243
           Iat=Iat+1
1244
        END DO
1245
     END DO
1246
     DEALLOCATE(NbAtType)
1247

  
1248
     NbTypes=NbType
1249

  
1250
     READ(IOIN,'(A)') Vasp_comment
1251
     READ(IOIN,'(A)') Vasp_direct
1252
     V_direct=Adjustl(Vasp_direct)
1253
     Call UpCase(V_direct)
1254
! PFL 2011 Mar 8 ->
1255
! We have to read the FFF flags :
1256
     DO I=1,Nat
1257
        READ(IOIN,*) Xtmp,ytmp,ztmp,FFF(1:3,I)
1258
        DO J=1,3
1259
           FFF(J,I)=AdjustL(FFF(J,I))
1260
           CALL Upcase(FFF(J,I))
1261
        END DO
1262
     END DO
1263
! <- PFL 2011 Mar 8 
1264

  
1265
  END IF
1266

  
1267
  ! In the case of VASP there is always the problem  of moving from one side
1268
  ! of the box to the other...
1269
! PFL 2011 Mar 14 ->
1270
! In fact, we have to make this check as soon as either
1271
! Input and/or Prog = VASP
1272
! PFL 2010 Dec 2 ->
1273
! Here we deal with input and not prog in fact 
1274
!  if (prog.EQ.'VASP') THEN
1275
  if ((Input.EQ.'VASP').OR.(Prog.EQ.'VASP')) THEN
1276
! <- PFL 2010 Dec 2
1277
! <- PFL 2011 Mar 14
1278
     Renum=.TRUE.
1279

  
1280
     ! V_direct has been set in Read_geom
1281
     IF (V_direct(1:6).EQ.'DIRECT') THEN
1282
        Latr(1:3,1)=Lat_a
1283
        Latr(1:3,2)=Lat_b
1284
        Latr(1:3,3)=Lat_c
1285
        B=1.
1286
        CALL Gaussj(Latr,3,3,B,1,1)
1287
     ELSE
1288
        Latr=0.
1289
        Latr(1,1)=1.d0
1290
        Latr(2,2)=1.d0
1291
        Latr(3,3)=1.d0
1292
     END IF
1293

  
1294
     ! Actualization of Frozen using the FFFF... 
1295
     ! Frozen has the advantage ie if given, it imposes _ALL_ the FFF flags.
1296
     IF (Frozen(1).NE.0) THEN
1297
        WRITE(IOOUT,*) "Frozen is given. Flags of the given POSCAR are overriden"
1298
        FFF='T'
1299

  
1300
        NFroz=0
1301
        DO WHILE (Frozen(NFroz+1).NE.0)
1302
           NFroz=NFroz+1
1303
           FFF(1:3,Frozen(NFroz))='F'
1304
        END DO
1305
     ELSE
1306
        WRITE(IOOUT,*) "Frozen not given : using  Flags of the given POSCAR"
1307
        NFroz=0
1308
        Frozen=0
1309
        DO I=1,Nat
1310
           IF ((FFF(1,I).EQ.'F').OR.(FFF(2,I).EQ.'F').OR.(FFF(3,I).EQ.'F')) THEN
1311
              FFF(1:3,I)='F'
1312
              NFroz=NFroz+1
1313
              Frozen(NFroz)=I
1314
           END IF
1315
        END DO
1316
        WRITE(IOOUT,*) "Frozen atoms:",Frozen(1:NFroz)
1317
     END IF
1318

  
1319

  
1320
  END IF
1321

  
1322
  IF (Prog=="VASP") THEN
1323
     IF (Vmd) THEN
1324
        if (debug) WRITE(*,*) "DBG MAIN, L803, VMD=T,NbTypes,AtTypes",NbTypes,AtTypes(1:NbTypes)
1325
        Line=""
1326
        DO I=1,NbTypes
1327
           Line=TRIM(Line) // " " // TRIM(AdjustL(AtTypes(I))) 
1328
        END DO
1329
        Vasp_Title=Trim(Line) // " " // Trim(adjustL(Vasp_Title))
1330
        if (debug) WRITE(*,*) "DBG MAIN, L809, VMD=T, Vasp_Title=",TRIM(Vasp_Title)
1331
     END IF
1332
     Call CheckPeriodicBound
1333
  END IF
1334

  
1335 902
  IF (COORD=="MIXED") Call TestCart(AutoCart)
1336 903

  
1337 904
  SELECT CASE(NFroz)
......
1450 1017
     END IF
1451 1018
  END IF
1452 1019

  
1453
  if (debug) THEN
1020
!  if (debug) THEN
1454 1021
     !Print *, 'L968, IntTangent(1,:)='
1455 1022
     !Print *, IntTangent(1,:)
1456
  END IF
1023
!  END IF
1457 1024

  
1458 1025
  ! We have read everything we needed... time to work :-)
1459 1026
  IOpt=0
......
1487 1054
!!! CAUTION : PFL 29.AUG.2008 ->
1488 1055
        ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
1489 1056
        ! This should be made  consistent !!!!!!!!!!!!!!!
1490
        GeomTmp=Reshape(reshape(XyzGeomF(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1057
        GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1491 1058
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1492 1059
!!! <- CAUTION : PFL 29.AUG.2008
1493 1060
     CASE ('ZMAT','MIXED')

Also available in: Unified diff