Révision 8 src/Path.f90

Path.f90 (revision 8)
331 331
     END subroutine Hinvup_BFGS_new
332 332

  
333 333

  
334
    FUNCTION InString(Line,String,Case,Clean,Back)  Result(Pos)
335

  
336
      Use VarTypes
337

  
338
      implicit none
339
! Input
340
      CHARACTER(*), INTENT(IN) :: Line
341
      CHARACTER(*), INTENT(IN) :: String
342
      LOGICAL, OPTIONAL, INTENT(IN) :: Case
343
      LOGICAL, OPTIONAL, INTENT(IN) :: Back
344
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean
345

  
346
! Output
347
! the position of String in Line (the first one) unless Back is present
348
      INTEGER(KINT) :: Pos
349
    END FUNCTION InString
350

  
351
    SUBROUTINE die(routine, msg, file, line, unit)
352

  
353
      Use VarTypes
354
      Use io_module
355

  
356
      implicit none
357
      character(len=*), intent(in)           :: routine, msg
358
      character(len=*), intent(in), optional :: file
359
      integer(KINT), intent(in), optional      :: line, unit
360

  
361
    END SUBROUTINE die
362

  
363
    SUBROUTINE Warning(routine, msg, file, line, unit)
364

  
365
      Use VarTypes
366
      Use io_module
367

  
368
      implicit none
369

  
370
      character(len=*), intent(in)           :: routine, msg
371
      character(len=*), intent(in), optional :: file
372
      integer(KINT), intent(in), optional      :: line, unit
373

  
374
    END SUBROUTINE Warning
375

  
376

  
334 377
  END INTERFACE
335 378

  
336 379
  Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, &
......
340 383
       OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, &
341 384
       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &
342 385
       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &
343
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile
386
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile,GplotFile
344 387

  
345 388
  Namelist/cartlist/list
346 389
  Namelist/frozenlist/list
......
355 398
  !
356 399
  WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 2007-2013"
357 400

  
358
  ! We read the files names
401
  ! We read the files  names
359 402
!  SELECT CASE (iargc())
360 403
  SELECT CASE (command_argument_count())
361 404
  CASE (2)
......
540 583
  GeomFile="EMPTY"
541 584
  AnaGeom=.FALSE.
542 585
  AnaFile="EMPTY"
586
  GplotFile="EMPTY"
543 587
  Nb=0
544 588
  NbCom=0
545 589
  PathOnly=.FALSE.
......
638 682
     if (CalcName=="EMPTY") CalcName="Path"
639 683
     if (ISuffix=="EMPTY")  ISuffix="com"
640 684
     if (OSuffix=="EMPTY")  OSuffix="log"
685
     UnitProg="au"
686
     ConvE=au2kcal
641 687
    CASE ("G03")
642 688
     Prog='GAUSSIAN'
643 689
     If (ProgExe=="EMPTY") ProgExe="g03"
644 690
     if (CalcName=="EMPTY") CalcName="Path"
645 691
     if (ISuffix=="EMPTY")  ISuffix="com"
646 692
     if (OSuffix=="EMPTY")  OSuffix="log"
693
     UnitProg="au"
694
     ConvE=au2kcal
647 695
    CASE ("EXT")
648 696
     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"
649 697
     if (CalcName=="EMPTY") CalcName="Ext"
650 698
     if (ISuffix=="EMPTY")  ISuffix="in"
651 699
     if (OSuffix=="EMPTY")  OSuffix="out"
700
     UnitProg="au"
701
     ConvE=au2kcal
652 702
    CASE ('MOPAC')
653 703
     If (ProgExe=="EMPTY") ProgExe="mopac"
654 704
     if (CalcName=="EMPTY") CalcName="Path"
655 705
     if (ISuffix=="EMPTY")  ISuffix="mop"
656 706
     if (OSuffix=="EMPTY")  OSuffix="out"
707
     UnitProg="au"
708
     ConvE=au2kcal
657 709
    CASE ("SIESTA")
658 710
     If (ProgExe=="EMPTY") ProgExe="siesta"
659 711
     if (CalcName=="EMPTY") CalcName="Path"
660 712
     if (ISuffix=="EMPTY")  ISuffix="fdf"
661 713
     if (OSuffix=="EMPTY")  OSuffix="out"
714
     UnitProg="eV"
715
     ConvE=eV2kcal
662 716
    CASE ("VASP")
663 717
     If (ProgExe=="EMPTY") ProgExe="VASP"
718
     UnitProg="eV"
719
     ConvE=eV2kcal
664 720
    CASE ("TM","TURBOMOLE")
665 721
     Prog="TURBOMOLE"
666 722
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
723
     UnitProg="au"
724
     ConvE=au2kcal
667 725
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
668 726
     ProgExe=""
727
     UnitProg="au"
728
     ConvE=au2kcal
669 729
    CASE DEFAULT
670 730
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
671 731
     STOP
......
750 810
     READ(IOIN,AnaList,IOSTAT=IoS)
751 811
     IF (IOS==0) THEN ! we have a AnaList namelist 
752 812
        Call ReadAnaList
813
     ELSE
814
! No AnaList namelist, we will print only the energy
815
        FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)'
753 816
     END IF
754 817
     IF (AnaFile=="EMPTY") THEN
755 818
        AnaFile=TRIM(PathName) // '.dat'
756 819
     END IF
757 820
     OPEN(IODat,File=AnaFile)
821
     IF (GplotFile=="EMPTY") THEN
822
        GplotFile=TRIM(PathName) // '.gplot'
823
     END IF
824
     
825
     OPEN(IODat,File=AnaFile)
826
     If (OptGeom<=0) THEN
827
        OPEN(IOGplot,File=GPlotFile)
828
        WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot -persist"
829
        WRITE(IoGplot,'(A)') " set pointsize 2"
830
        WRITE(IoGplot,'(A)') " xcol=1"
831
        WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3
832
     END IF
758 833
  END IF
759 834

  
760 835
        
......
799 874
  AtName=""
800 875
  MassAt=1.
801 876

  
877
! As we have used rewind many times, the position in the input file
878
! is blurry (at least !) we thus have to go back to the Geometries description
879
! (TO DO: use a better Input parser !)
880

  
881
  REWIND(IOIN)
882
! We search the '&path' namelist
883
  READ(IOIN,'(A)',Iostat=Ios) Line
884
  Fini=.FALSE.
885
  DO WHILE ((Ios==0).AND.InString(Line,"&PATH")==0)
886
     if (debug) WRITE(*,*) "Lookgin for Path:",TRIM(Line)
887
     READ(IOIN,'(A)',Iostat=Ios) Line
888
  END DO
889
  if (debug) WRITE(*,*) "Path found:",TRIM(LINE)
890
  ! We are on the &path line, we move to the end
891
  Call NoComment(Line)
892
  DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
893
     if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
894
     READ(IOIN,'(A)',Iostat=Ios) Line
895
     Call NoComment(Line)
896
     Call noString(Line)
897
  END DO
898
! We are at then end of &Path / namelist
899
! We scan for other Namelists
900
  READ(IOIN,'(A)',Iostat=Ios) Line
901
  if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
902
  DO WHILE  ((IoS==0).AND.(Index(Line,'&')/=0))
903
     Call NoComment(Line)
904
     if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
905
     Idx=InString(Line,'Analist')
906
     DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
907
        if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
908
        READ(IOIN,'(A)',Iostat=Ios) Line
909
        Call NoComment(Line)
910
     END DO
911
     If (Idx/=0) THEN
912
! we have just read &Analist namelist, we have to skip the variables description
913
        DO I=1,Nb
914
           READ(IOIN,'(A)',Iostat=Ios) Line
915
           if (debug) WRITe(*,*) "We skip Analist",TRIM(LINE)
916
        END DO
917
     END IF
918
     READ(IOIN,'(A)',Iostat=Ios) Line
919
     if (debug) WRITe(*,*) "We skip the line:",TRIM(LINE)
920
  END DO
921
  BACKSPACE(IOIN)
922

  
802 923
  ! We read the initial geometries
803 924
  Call Read_Geom(input)
804 925

  
......
809 930
  END DO
810 931

  
811 932
 If (AnaGeom) THEN
933
    If (NbVar>0) THEN
812 934
! We print the description of the varialbes in AnaFile
813
    Call PrintAnaList(IoDat)
935
       Call PrintAnaList(IoDat)
936
       If (OptGeom<=0) THEN
937
          WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
938
       ELSE
939
          WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
940
       END IF
941
    ELSE
942
       If (OptGeom<=0) THEN
943
          WRITE(IoDat,'(A)') "#  s    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"
944
       ELSE
945
          WRITE(IoDat,'(A)') "#  Iter    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"       
946
       END IF
947
    END IF
814 948
 END IF
815 949

  
816 950
  ! PFL 23 Sep 2008 ->
......
879 1013
              Dist=Norm1
880 1014
           END IF
881 1015
        END DO
882
        IAt2=Frozen(2)
883 1016
        ListAtFroz(Iat2)=.TRUE.
884 1017
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2
885 1018
        x2(1)=x1(Iat2)-x1(Iat1)
......
888 1021
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
889 1022
        Dist=-1.
890 1023
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
1024
        Iat3=1
1025
        DO WHILE (ListAtFroz(Iat3).AND.(Iat3<=Nat))
1026
           Iat3=Iat3+1
1027
        END DO
891 1028
        DO I=1,Nat
892 1029
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
893 1030
           x2(2)=x1(Iat2)-x1(I)
......
909 1046
     DO I=2,NGeomI
910 1047
        ! First, we align the Ith geometry on the first one, using only
911 1048
        ! the positions of the "frozen" atoms. (see above: it might be that
912
        ! ListAtFroz contains non frozen atoms usd to align the geometries
1049
        ! ListAtFroz contains non frozen atoms used to align the geometries
913 1050
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
914 1051
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
915 1052
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
916 1053
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
917 1054

  
918

  
919 1055
        FChkFrozen=.FALSE.
920 1056
        DO J=1,NFroz
921 1057
           IAt=Frozen(J)
......
2101 2237
  Close(IOIN)
2102 2238
  Close(IOOUT)
2103 2239
  IF (AnaGeom) Close(IODAT)
2240
  IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot)
2104 2241

  
2105 2242
END PROGRAM PathOpt

Formats disponibles : Unified diff