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 MayJune > 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 20072011"


342 
WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 20072013"


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 MurtaghSargent" 
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=N3at6 ! 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 userdefined 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