-
Notifications
You must be signed in to change notification settings - Fork 0
/
BML110.FORT10
1929 lines (1929 loc) · 51.5 KB
/
BML110.FORT10
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
C**********************************************************************
C Library for BMATIN6 fortran.
C**********************************************************************
c Moved to PSI area on 2/3/89 by clj.
C**********************************************************************
C NOTIFICATION OF PROGRAM CHANGE: *
C BY: TRACY HAMILTON *
C DATE: JULY 18,1988 *
C REASON: SCALE FORCES CORRECTLY IN DO 20 LOOP *
C**********************************************************************
SUBROUTINE EFC(NUMIT,NQ,F1,QQ,FI,C,MODEBM,QQ1,IFSTRE,GEC)
IMPLICIT REAL*8 (A-H,O-Z)
C
C THIS ROUTINE SETS UP ALL THE COMMON BLOCKS NEEDED PRIOR TO CALLING
C LINK 110.
C
COMMON/MODOUT/VMODE1(50)
COMMON/IOP/IOP(50)
COMMON/IO/IN,IOUT,INP2,IPUN,IGMUP
COMMON/GRDNT/ENERGY,F(50),FRCNST(1275),NVAR,IGETFC
COMMON/OPTEF/X(50),XNAME(50),OLDF(50),D(50),HESS(50,50),
$ IC(50),VMODE(50),MODE,ES,NSTEP,MXSTEP,IHESS,IS,
$ NEGREQ,DMAX,DD,CONVF,CNVFMX,CONVX,CNVXMX,IUPDAT,
$ MXHESS,EIGMAX,EIGMIN,PRNT,NRFLAG,DUMP,IDUM
COMMON/ASTRE/STRE(50)
DIMENSION F1(1),QQ(1),FI(1),C(1),QQ1(1),CORE(1),IFSTRE(1),GEC(1)
LOGICAL STRE,IFSTRE
DATA TOANG/0.52917706D0/, HARTRE/4.359813653D0/
CONV1 = TOANG/HARTRE
CONV2 = TOANG*CONV1
DO 5 I=1,NQ
5 STRE(I) = IFSTRE(I)
C
NSTEP = NUMIT
MODE = MODEBM
C
C NUMBER OF SYMMETRY INTERNAL COORDINATES:
NVAR = NQ
C
C VALUES OF SYMMETRY INTERNAL COORDINATES:
DO 10 I=1,NVAR
X(I) = QQ(I)
C CONVERT ANG --> BOHR
IF(STRE(I)) X(I) = X(I) / TOANG
10 CONTINUE
C WRITE(6,666) (X(I),I=1,NVAR)
C 666 FORMAT('X ',6F10.7)
C
C FORCES (CONVERT MDYNE --> HARTREE/BOHR):
DO 20 I=1,NVAR
IF(STRE(I)) THEN
F(I) = FI(I) * CONV1
ELSE
F(I) = FI(I) / HARTRE
ENDIF
20 CONTINUE
C WRITE(6,667) (F(I),I=1,NVAR)
C 667 FORMAT('F ',6F10.7)
C
C HESSIAN:
IJ = 0
DO 31 I=1,NVAR
DO 30 J=1,NVAR
IJ = IJ+1
HESS(I,J) = C(IJ)
IF(STRE(I).AND.STRE(J)) THEN
C MDYNE/ANG --> HARTREE/BOHR**2
HESS(I,J) = HESS(I,J) * CONV2
ELSE IF(STRE(I).OR.STRE(J)) THEN
C MDYNE/RAD --> HARTREE/BOHR/RAD
HESS(I,J) = HESS(I,J) * CONV1
ELSE
C MDYNE*ANG/RAD**2 --> HARTREE/RAD**2
HESS(I,J) = HESS(I,J) / HARTRE
ENDIF
30 CONTINUE
C WRITE(6,668) (HESS(I,J),J=1,NVAR)
C 668 FORMAT('HESS ',6F10.7)
31 CONTINUE
C
IF(NSTEP.GT.0) THEN
C
C OLD GRADIENT (NEGATIVE OF OLD FORCES) (CONVERT TO HARTREE/BOHR):
DO 40 I=1,NVAR
40 OLDF(I) = -F1(I) * CONV1
C WRITE(6,669) (OLDF(I),I=1,NVAR)
C 669 FORMAT('OLDF ',6F10.7)
C
C SKIP OVER FIRST LINE IN INP4:
CTPH READ(INP4,1000) DD
C
C OLD STEP SIZE AND STEP:
CTPH READ(INP4,1000) DD
CTPH READ(INP4,1000) (D(I),I=1,NVAR)
DD = 0.0D0
DO 45 I=1,NVAR
D(I) = GEC(I)
IF(STRE(I)) D(I) = D(I) / TOANG
DD = DD + D(I)**2
45 CONTINUE
DD = DSQRT(DD)
C WRITE(6,670) (D(I),I=1,NVAR)
C 670 FORMAT('D ',6F10.7)
C
C POSSIBLY RECOVER OLD HESSIAN MODE:
IF(MODE.GT.0) THEN
C READ(INP3,1010) MODE
C READ(INP3,1020) (VMODE(I),I=1,NVAR)
DO 46 I=1,NVAR
READ(INP2,1020) VMODE(I)
46 CONTINUE
IF(MODE.EQ.1) MODE=0
ENDIF
ENDIF
C
C
C
CALL OPTEFC(CORE)
C
C
C
CTPH REWIND INP4
CTPH WRITE(INP4,1010) -1,NSTEP
C
C SAVE STEP SIZE AND STEP:
CTPH WRITE(INP4,1000) DD
CTPH WRITE(INP4,1000) (D(I),I=1,NVAR)
C
C SAVE MODE:
CTPH REWIND INP3
CTPH WRITE(INP3,1010) MODE
CTPH WRITE(INP3,1020) (VMODE(I),I=1,NVAR)
C
C PUT HESSIAN BACK INTO BMAT FORM:
MODEBM = MODE
DO 51 I=1,NVAR
VMODE1(I) = VMODE(I)
51 CONTINUE
IJ = 0
DO 70 I=1,NVAR
DO 70 J=1,NVAR
IJ = IJ+1
IF(STRE(I).AND.STRE(J)) THEN
C HARTREE/BOHR**2 --> MDYNE/ANG
HESS(I,J) = HESS(I,J) / CONV2
ELSE IF(STRE(I).OR.STRE(J)) THEN
C HARTREE/BOHR/RAD --> MDYNE/RAD
HESS(I,J) = HESS(I,J) / CONV1
ELSE
C HARTREE/RAD**2 --> MDYNE*ANG/RAD**2
HESS(I,J) = HESS(I,J) * HARTRE
ENDIF
70 C(IJ) = HESS(I,J)
C WRITE(6,671) (C(I),I=1,NVAR*NVAR)
C 671 FORMAT('C ',6F10.7)
C
C COPY DISPLACEMENTS INTO BMAT ARRAY (AND CONVERT BOHR --> ANG):
DO 80 I=1,NVAR
QQ1(I) = D(I)
IF(STRE(I)) QQ1(I) = QQ1(I) * TOANG
80 CONTINUE
C
RETURN
C
1000 FORMAT(F10.6)
1010 FORMAT(2I5)
1020 FORMAT(8F10.6)
C
END
SUBROUTINE OPTEFC(CORE)
IMPLICIT REAL*8(A-H,O-Z)
C
C
C **************************************************************
C * *
C * DRIVER FOR EIGENVECTOR FOLLOWING TRANSITION STATE SEARCH *
C * BASED ON "ON FINDING TRANSITION STATES" BY CERJAN AND *
C * MILLER (J.CHEM.PHYS. 75 (1981) 2800); "WALKING ON *
C * POTENTIAL ENERGY SURFACES" BY SIMONS, JORGENSON, TAYLOR *
C * AND OZMENT (J.PHYS.CHEM. 87 (1983) 2745) AND "SEARCH *
C * FOR STATIONARY POINTS ON SURFACES" BY BANERJEE, ADAMS, *
C * SIMONS AND SHEPARD (J.PHYS.CHEM. 89 (1985) 52). *
C * *
C * ORIGINAL VAX CODE BY JON BAKER JB MARCH 1985 *
C * *
C **************************************************************
C
C
C OPTIONS COMMON/IOP/
C
C IOP(5) NATURE OF REQUIRED STATIONARY POINT
C 0 FIND A TS (DEFAULT)
C 1 FIND A MINIMUM
C
C IOP(6) MAXIMUM NUMBER OF STEPS ALLOWED
C 0 MIN(40,NVAR+20) (DEFAULT)
C (WHERE NVAR = NUMBER OF VARIABLES)
C N N STEPS
C
C IOP(7) CONVERGENCE CRITERIA ON RMS GRADIENT
C 0 0.0003 (DEFAULT)
C N 0.001/N
C NOTE: THE OTHER CONVERGENCE CRITERIA ARE
C MAXIMUM GRADIENT = 1.5 * RMS GRADIENT
C RMS DISPLACEMENT = 4.0 * RMS GRADIENT
C MAX DISPLACEMENT = 1.5 * RMS DISPLACEMENT
C
C IOP(8) MAXIMUM STEPSIZE ALLOWED DURING OPTIMIZATION
C 0 DMAX = 0.3 (DEFAULT)
C N DMAX = 0.01*N
C
C IOP(10) INPUT OF INITIAL SECOND DERIVATIVE MATRIX
C ALL VALUES MUST BE IN ATOMIC UNITS
C 0 ESTIMATE ENTIRE HESSIAN NUMERICALLY (TS DEFAULT)
C ESTIMATE DIAGONAL ELEMENTS ONLY (MIN DEFAULT)
C 1 READ IN FULL FRCNST MATRIX (LOWER TRIANGLE)
C FORMAT (8F10.6)
C 2 READ IN SELECTED ELEMENTS OF FRCNST
C READ I,J,FRCNST(I,J) 2I3, F20.0
C 3 READ FRCNST MATRIX FROM THE CHECKPOINT FILE
C 4 CALCULATE HESSIAN ANALYTICALLY
C 5 READ CARTESIAN FORCE CONSTANTS FROM CHECKPOINT FILE
C
C IOP(11) HESSIAN RECALCULATION
C 0 UPDATE THE HESSIAN ONLY (DEFAULT)
C N RECALCULATE THE HESSIAN (ANALYTICALLY IF
C POSSIBLE, ELSE NUMERICALLY) EVERY N POINTS
C (OBS! IF THIS OPTION IS INVOKED, RESTARTS FROM THE
C CHECKPOINT FILE SHOULD NOT BE ATTEMPTED EXCEPT VIA A
C NON-STANDARD ROUTE. EXCEPTION: ANALYTICAL DERIVATIVES
C AT ALL POINTS VIA "CALCALL" CAN BE RESTARTED)
C
C IOP(13) TYPE OF HESSIAN UPDATE
C 0 POWELL UPDATE (DEFAULT)
C 1 BFGS UPDATE (USED IN CONJUNCTION WITH MIN)
C 2 BFGS UPDATE WITH SAFEGUARDS TO ENSURE RETENTION
C OF POSITIVE DEFINITENESS
C
C IOP(16) MAXIMUM ALLOWABLE MAGNITUDE OF HESSIAN EIGENVALUES
C IF THIS MAGNITUDE IS EXCEEDED, THE EIGENVALUE IS REPLACED
C 0 EIGMAX = 25.0 (DEFAULT)
C N EIGMAX = 0.1*N
C
C IOP(17) MINIMUM ALLOWABLE MAGNITUDE OF HESSIAN EIGENVALUES
C SIMILAR TO IOP(16)
C 0 EIGMIN = 0.0001 (DEFAULT)
C N EIGMIN = 1.0/N
C
C IOP(19) SEARCH SELECTION
C 0 P-RFO OR RFO STEP ONLY (DEFAULT)
C 1 P-RFO OR RFO STEP FOR "WRONG" HESSIAN
C OTHERWISE NEWTON-RAPHSON
C
C IOP(21) EXPERT SWITCH
C 0 NORMAL MODE (DEFAULT)
C 1 EXPERT MODE
C CERTAIN CUTOFFS TO CONTROL OPTIMIZATION RELAXED
C
C IOP(33) PRINT OPTION
C 0 ON (DEFAULT)
C 1 OFF TURNS OFF EXTRA PRINTING
C (DEFAULT OF "ON" BY POPULAR REQUEST)
C
C IOP(34) DUMP OPTION
C 0 OFF (DEFAULT)
C 1 ON TURNS ON DEBUG PRINTING
C
C IOP(35) RESTART OPTION
C 0 NORMAL OPTIMIZATION (DEFAULT)
C 1 FIRST POINT OF A RESTART
C RECOVER GEOMETRY ETC.. FROM CHECKPOINT FILE
C
C IOP(36) CHECKPOINT OPTION
C 0 CHECKPOINTING AS NORMAL (DEFAULT)
C 1 SUPPRESS CHECKPOINTING
C
C IOP(37) D2E CLEANUP OPTION
C 0 NO SPECIAL ACTION TAKEN (DEFAULT)
C 1 THIS IS THE LAST POINT AT WHICH ANALYTICAL
C DERIVATIVES ARE AVAILABLE; CLEAN UP THE RWF
C
C
C MODE FOLLOWING: MODE FOLLOWING IS TURNED ON VIA THE IC ARRAY, WHICH
C IS INPUT WITH THE Z-MATRIX VARIABLES IN LINK 101. ADDING 4 TO THE IC
C ENTRY FOR A PARTICULAR VARIABLE WILL CAUSE A TRANSITION STATE SEARCH
C TO FOLLOW THE HESSIAN MODE WITH THE LARGEST MAGNITUDE COMPONENT FOR
C THAT VARIABLE. ADDING 10 TO THE IC ENTRY FOR THE KTH VARIABLE WILL
C FOLLOW THE KTH HESSIAN MODE. ONLY ONE IC ENTRY SHOULD BE MODIFIED IN
C THIS WAY I.E. ONLY ONE MODE SHOULD BE FOLLOWED AT A TIME.
C
C
c REAL*8 CORE(1)
REAL*8 U(50,50),HESSC(50,50),EIGVAL(50),FX(50),TVEC(50),tvec2(50)
LOGICAL PRNT,NRFLAG,DUMP,CNVGRD,last
COMMON/IOP/IOP(50)
COMMON/IO/IN,IOUT
COMMON/GRDNT/ENERGY,F(50),FRCNST(1275),NVAR,IGETFC
COMMON/OPTEF/X(50),XNAME(50),OLDF(50),D(50),HESS(50,50),
$ IC(50),VMODE(50),MODE,ES,NSTEP,MXSTEP,IHESS,IS,
$ NEGREQ,DMAX,DD,CONVF,CNVFMX,CONVX,CNVXMX,IUPDAT,
$ MXHESS,EIGMAX,EIGMIN,PRNT,NRFLAG,DUMP,IDUM
c COMMON/ZSUBST/ANAMES(50),VALUES(50),INTVEC(50),FPVEC(50)
EQUIVALENCE (FX(1),TVEC(1))
C
c DATA IGRDNT/511/,LGRDNT/1327/, IZSUBS/570/,LZSUBS/175/
c DATA IOPTEF/575/,LOPTEF/2790/
DATA ZERO/0.0D0/, ONE/1.0D0/
C
C
C INITIALIZE THIS LINK
C THE ILSW BIT WILL BE SET UNLESS THIS IS THE FIRST ENTRY TO
C LINK 110, IN WHICH CASE WE'LL JUST DO SOME INITIALIZATION
C
c CALL DRUM
WRITE(IOUT,1000)
c CALL ILSW(2,23,IGRD)
c CALL ILSW(1,23,1)
c IF(IGRD.EQ.1.OR.IOP(35).GT.0) GO TO 100
C
C READ IN DATA AND INITIALIZE OPTIMIZATION
C
IF(IOP(5).EQ.0) WRITE(IOUT,1005)
IF(IOP(5).NE.0) WRITE(IOUT,1006)
c WRITE(IOUT,1010)
CALL INITEF(CORE)
c CALL EXITEF(0)
c RETURN
C
C START LOOP FOR SECOND AND LATER CALLS TO LINK 110
C
C RESTORE THE CHECKPOINTED INFORMATION IF THIS IS THE FIRST
C POINT OF A RESTART FROM THE CHECKPOINT FILE
C
c 100 IF(IOP(35).EQ.1) CALL CHKPNT(1)
C
C RECOVER COMMON BLOCKS FROM RWF
C
c CALL TREAD(IZSUBS,ANAMES,LZSUBS,1,LZSUBS,1,0)
c CALL TREAD(IGRDNT,ENERGY,LGRDNT,1,LGRDNT,1,0)
c CALL TREAD(IOPTEF,X,LOPTEF,1,LOPTEF,1,0)
c IF(NEGREQ.EQ.1) WRITE(IOUT,1005)
c IF(NEGREQ.EQ.0) WRITE(IOUT,1006)
c IF(IOP(6).NE.0) MXSTEP=IOP(6)
c IF(IOP(36).EQ.0) CALL CHKPNT(0)
C
C CONVERT FORCES TO GRADIENT
C
CALL DNEGV(F,NVAR)
C
C CHECK ON STATUS OF HESSIAN
C THIS IS GIVEN BY THE INTEGER FLAG IHESS
C IHESS=0 HESSIAN UNAVAILABLE
C (BEING ESTIMATED NUMERICALLY)
C IHESS=1 FIRST POINT WITH CONSTRUCTED HESSIAN
C IHESS=N NTH POINT WITH HESSIAN
C (SHOULD BE UPDATED PRIOR TO USE)
C
c IF(IHESS.LE.0) CALL HESSEF
C
C IF HESSIAN STILL UNAVAILABLE UPDATE PARAMETER LIST AND EXIT
C
c IF(IHESS.EQ.0) THEN
c WRITE(IOUT,1015)
c WRITE(IOUT,1020) NSTEP,IS
c CALL EXITEF(0)
c RETURN
c ENDIF
C
C UPDATE THE HESSIAN (UNLESS FIRST POINT)
c (or unless already updated by BMAT or calculated analytically)
C
c IF(IHESS.GT.1) CALL UPDATE(EIGVAL,TVEC)
IF(nstep.gt.0 .and. iop(13).ge.0) CALL UPDATE(EIGVAL,TVEC)
C
C IF THIS IS THE LAST POINT AT WHICH ANALYTICAL SECOND
C DERIVATIVES WILL BE AVAILABLE, CLEAN UP THE RWF
C
c IF(IOP(37).NE.0) CALL CLND2E
C
C HESSIAN READY FOR USE
C START OPTIMIZATION PROPER
C
NSTEP=NSTEP+1
WRITE(IOUT,1025) NSTEP
c IHESS=IHESS+1
C
C DIAGONALIZE HESSIAN
C FIRST COPY HESS INTO HESSC SINCE HESSIAN IS DESTROYED
C BY DIAGONALIZATION ROUTINE
C
c CALL DMCOPY(HESSC,HESS,NVAR)
c CALL DODIAG(NVAR,50,HESSC,U,EIGVAL,TVEC)
call cntrct(hessc,hess,nvar)
call rsp(50,nvar,2500,hessc,eigval,1,u,tvec,tvec2)
C
C PRINT THE HESSIAN, ITS EIGENVECTORS AND EIGENVALUES IF
C REQUESTED. VECTORS ARE IN U, EIGENVALUES IN EIGVAL
C
IF(PRNT) THEN
IF(DUMP) WRITE(IOUT,1030)
IF(DUMP) CALL DMPRNT(HESS,NVAR,IOUT)
WRITE(IOUT,1035)
c CALL MATPRT(U,50,50,NVAR,NVAR,2,0,XNAME,XNAME,0,EIGVAL,1)
last = .false.
nstar = 1
10 nend = nstar + 5
if(nend.ge.nvar) then
last = .true.
nend = nvar
endif
write(iout,1045) (eigval(l),l=nstar,nend)
write(iout,1045)
do 997 i=1,nvar
write(iout,1045) (u(i,l),l=nstar,nend)
997 continue
write(iout,1045)
nstar = nend + 1
if(.not.last) go to 10
C
ELSE
C
WRITE(IOUT,1040)
WRITE(IOUT,1045) (EIGVAL(L),L=1,NVAR)
ENDIF
C
C CHECK FOR MAXIMUM AND MINIMUM ALLOWED MAGNITUDES
C FOR THE HESSIAN EIGENVALUES
C
CALL MAGCHK(EIGVAL,NVAR)
C
C CALCULATE NEG, THE NUMBER OF NEGATIVE EIGENVALUES
C
CALL FNDNEG(EIGVAL,NEG,NVAR)
C
C FORM THE FX VECTOR
C (THE COMPONENT OF F ALONG THE LOCAL EIGENDIRECTIONS U)
C
DO 20 I=1,NVAR
FX(I)=DVTV(U(1,I),F,NVAR)
20 CONTINUE
C
C IF WE ARE IN THE "RIGHT" REGION OF THE ENERGY SURFACE
C I.E. THE HESSIAN HAS THE REQUIRED NUMBER OF NEGATIVE EIGENVALUES
C AND NRFLAG IS ON, TAKE THE NEWTON-RAPHSON STEP
C
IF(NRFLAG.AND.NEG.EQ.NEGREQ) THEN
WRITE(IOUT,1050)
CALL FORMNR(U,EIGVAL,FX,D,NVAR)
C
ELSE
C
C TAKE THE P-RFO STEP FOR A TS SEARCH AND
C THE SIMPLE RFO STEP FOR A MINIMUM SEARCH
C
IF(NEG.EQ.NEGREQ.AND.NEGREQ.EQ.1) WRITE(IOUT,1055)
IF(NEG.EQ.NEGREQ.AND.NEGREQ.EQ.0) WRITE(IOUT,1056)
IF(NEG.NE.NEGREQ.AND.NEGREQ.EQ.1) WRITE(IOUT,1057)
IF(NEG.NE.NEGREQ.AND.NEGREQ.EQ.0) WRITE(IOUT,1058)
CALL FORMD(U,EIGVAL,FX,NVAR)
ENDIF
C
C WE NOW HAVE A NEW STEP IN D
C CHECK THAT THE STEPSIZE DOES NOT EXCEED DMAX
C IF SO, SCALE
C
CALL CHECKD(D,DMAX,DD,SKAL,NVAR)
C
C PREDICT THE CHANGE IN ENERGY
C
CALL ESTIME(EIGVAL,FX,DD,SKAL,CHNGE,NVAR)
C
C TEST FOR CONVERGENCE AND PRINT OUT CURRENT PARAMETER VALUES
C
CALL CONVEF(NEG,CNVGRD,TVEC)
WRITE(IOUT,1060) CHNGE
C
C IF CONVERGED THEN EXIT
C
IF(CNVGRD) THEN
c CALL EXITEF(1)
write(iout,1000)
RETURN
ENDIF
C
C FORM THE NEW X VECTOR AND SAVE THE CURRENT GRADIENT
C
c CALL DVADDV(X,X,D,NVAR)
c CALL DVCOPY(OLDF,F,NVAR)
C
C CHECK IF THE HESSIAN IS TO BE RECALCULATED/REESTIMATED
C
c IF(IHESS.GT.MXHESS) THEN
c IHESS=0
c IS=0
c CALL GET202(ICHAIN,1)
c CALL EXITEF(ICHAIN)
c RETURN
c ENDIF
C
C AND EXIT
C
c CALL EXITEF(0)
write(iout,1000)
C
RETURN
C
1000 FORMAT(/24('EF-'))
C
1005 FORMAT(//'EIGENVECTOR FOLLOWING TRANSITION STATE SEARCH')
C
1006 FORMAT(//'EIGENVECTOR FOLLOWING MINIMUM SEARCH')
C
1010 FORMAT(' INITIALIZATION PASS')
C
1015 FORMAT(//' NUMERICALLY ESTIMATING SECOND DERIVATIVES')
C
1020 FORMAT(' ITERATION ',I3,' VARIABLE ',I2)
C
1025 FORMAT(//'ITERATION ',I3)
C
1030 FORMAT(//' HESSIAN MATRIX')
C
1035 FORMAT('EIGENVECTORS AND EIGENVALUES'/)
C
1040 FORMAT(//'EIGENVALUES OF THE HESSIAN')
C
c1045 FORMAT(5X,6F12.6)
1045 FORMAT(6F12.6)
C
1050 FORMAT(' HESSIAN HAS REQUIRED LOCAL STRUCTURE'/
$ ' TAKING NEWTON-RAPHSON STEP')
C
1055 FORMAT(' TS SEARCH. TAKING P-RFO STEP')
C
1056 FORMAT(' MINIMUM SEARCH. TAKING SIMPLE RFO STEP')
C
1057 FORMAT(' HESSIAN DOES NOT HAVE THE DESIRED LOCAL STRUCTURE'/
$ ' TAKING P-RFO STEP')
C
1058 FORMAT(' HESSIAN DOES NOT HAVE THE DESIRED LOCAL STRUCTURE'/
$ ' TAKING SIMPLE RFO STEP')
C
1060 FORMAT(' PREDICTED CHANGE IN ENERGY ',F9.6)
C
END
SUBROUTINE CHECKD(D,DMAX,DD,SKAL,NVAR)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 D(1)
COMMON/IO/IN,IOUT
C
C GET CURRENT STEPSIZE
C
DD=DVTV(D,D,NVAR)
DD=DSQRT(DD)
C
IF(DD.GT.DMAX) THEN
SKAL=DMAX/DD
WRITE(IOUT,1000) SKAL
CALL DVSKAL(D,D,SKAL,NVAR)
DD=DMAX
ENDIF
C
WRITE(IOUT,1010) DD
RETURN
C
1000 FORMAT(' CALCULATED STEP TOO LARGE. STEP SCALED BY ',F9.6)
C
1010 FORMAT(' STEP TAKEN. STEPSIZE IS ',F9.6)
C
END
SUBROUTINE CONVEF(NEG,CNVGRD,TVEC)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION TVEC(1)
c DIMENSION TVEC(1),IRESLT(3)
character*3 ireslt
LOGICAL PRNT,NRFLAG,DUMP,CNVGRD,stre
common/astre/stre(50)
COMMON/IO/IN,IOUT
c COMMON/PHYCON/TOANG,PHYCON(29)
COMMON/GRDNT/ENERGY,F(50),FRCNST(1275),NVAR,IGETFC
COMMON/OPTEF/X(50),XNAME(50),OLDF(50),D(50),HESS(50,50),
$ IC(50),VMODE(50),MODE,ES,NSTEP,MXSTEP,IHESS,IS,
$ NEGREQ,DMAX,DD,CONVF,CNVFMX,CONVX,CNVXMX,IUPDAT,
$ MXHESS,EIGMAX,EIGMIN,PRNT,NRFLAG,DUMP,IDUM
c COMMON/ZMAT/IANZ(50),IZ(50,4),BL(50),ALPHA(50),BETA(50),
c $ LBL(50),LALPHA(50),LBETA(50),NZ,NVARRD
C
c DATA IOZMAT/507/,LZMAT/351/
data toang/0.52917706d0/, todeg/57.296d0/
C
C PRINT OUT CURRENT X, CURRENT GRADIENTS, DISPLACEMENTS
C AND NEWX
C
WRITE(IOUT,1000)
WRITE(IOUT,1010)
DO 20 I=1,NVAR
TVEC(I)=X(I)+D(I)
if(stre(i)) then
temp = tvec(i) * toang
else
temp = tvec(i) * todeg
endif
WRITE(IOUT,1020) I,X(I),F(I),D(I),TVEC(I),temp
20 CONTINUE
C
C FIND MAXIMUM AND RMS GRADIENTS AND DISPLACEMENTS
C
FMAX=DABS(F(1))
DXMAX=DABS(D(1))
IF(NVAR.GT.1) THEN
DO 30 I=2,NVAR
IF(DABS(F(I)).GT.FMAX) FMAX=DABS(F(I))
IF(DABS(D(I)).GT.DXMAX) DXMAX=DABS(D(I))
30 CONTINUE
ENDIF
C
FSQ=DVTV(F,F,NVAR)
FRMS=DSQRT(FSQ/NVAR)
DXRMS=DSQRT(DD*DD/NVAR)
C
C CHECK CONVERGENCE
C
WRITE(IOUT,1030)
CALL CONVGD(FMAX,CNVFMX,IRESLT)
WRITE(IOUT,1032) FMAX,CNVFMX,IRESLT
CALL CONVGD(FRMS,CONVF,IRESLT)
WRITE(IOUT,1034) FRMS,CONVF,IRESLT
CALL CONVGD(DXMAX,CNVXMX,IRESLT)
WRITE(IOUT,1036) DXMAX,CNVXMX,IRESLT
CALL CONVGD(DXRMS,CONVX,IRESLT)
WRITE(IOUT,1038) DXRMS,CONVX,IRESLT
C
C SET CONVERGENCE FLAG AND TEST FOR MAXIMUM ITERATIONS
C EXCEEDED
C
IF(FMAX.GT.CNVFMX.OR.FRMS.GT.CONVF.OR.
$ DXMAX.GT.CNVXMX.OR.DXRMS.GT.CONVX) THEN
CNVGRD=.FALSE.
c IF(NSTEP.GE.MXSTEP) THEN
c WRITE(IOUT,1040)
c CALL ILSW(1,25,0)
c CALL ILSW(1,27,1)
c CNVGRD=.TRUE.
c ENDIF
RETURN
ENDIF
C
C CONVERGED. PRINT FINAL SOLUTION
C
CNVGRD=.TRUE.
WRITE(IOUT,1050)
C
C TEST IF THE HESSIAN HAS THE REQUIRED NUMBER OF
C NEGATIVE EIGENVALUES
C
IF(NEG.NE.NEGREQ) WRITE(IOUT,1060)
C
C PRINT FINAL PARAMETERS
C
c CALL TREAD(IOZMAT,IANZ,LZMAT,1,LZMAT,1,0)
c DO 40 I=1,NVAR
c IC(I)=99
c 40 CONTINUE
C
c CALL PRMTBL(1,XNAME,X,IC,F,NVAR,LBL,NZ,TOANG)
C
RETURN
C
1000 FORMAT(//'CURRENT PARAMETER VALUES (INTERNAL COORDINATES)',
$ 17x,'Angstrom')
C
1010 FORMAT(3X,'I',8X,'X',9X,'GRADIENT',4X,
$ 'DISPLACEMENT',7X,'NEWX',5x,'or degrees')
C
1020 FORMAT(I4,4(2X,F10.6,2X),3x,f9.4)
C
1030 FORMAT(//8X,4HITEM,15X,5HVALUE,5X,9HTHRESHOLD,2X,
$ 10HCONVERGED?)
C
1032 FORMAT(/14H MAXIMUM FORCE,12X,F8.6,5X,F8.6,5X,a3)
C
1034 FORMAT(14H RMS FORCE,12X,F8.6,5X,F8.6,5X,a3)
C
1036 FORMAT(21H MAXIMUM DISPLACEMENT,5X,F8.6,5X,F8.6,5X,a3)
C
1038 FORMAT(21H RMS DISPLACEMENT,5X,F8.6,5X,F8.6,5X,a3)
C
1040 FORMAT(//' ***********************************'/
$ ' ** OPTIMIZATION STOPPED **'/
$ ' ** MAXIMUM ITERATIONS EXCEEDED **'/
$ ' ***********************************'//)
C
1050 FORMAT(//' *************************************************'/
$ ' ** CONVERGENCE CRITERIA APPARENTLY SATISFIED **'/
$ ' *************************************************'//)
C
1060 FORMAT(//'|| WARNING: THE HESSIAN HAS THE WRONG NUMBER OF',
$ ' NEGATIVE EIGENVALUES ||'/)
C
END
c SUBROUTINE ESTM
c IMPLICIT REAL*8(A-H,O-Z)
c LOGICAL PRNT,NRFLAG,DUMP
C
C MAKES GUESSES AT THE DIAGONAL SECOND DERIVATIVES
C THIS ROUTINE IS A MODIFIED VERSION OF THE GAUSSIAN 82 ROUTINE
C OF THE SAME NAME USED IN LINK 103.
C
C BENDING FORCE CONSTANTS ARE 1.35 FOR MINIMAL BASIS, ELSE 1.0.
C FOR STRETCHES THE VALUE CHOSEN DEPENDS UPON THE INTERNUCLEAR
C DISTANCE AND THE ROWS OF THE PERIODIC TABLE IN WHICH THE TWO
C AFFECTED ATOMS RESIDE AS WELL AS THE BASIS.
C DUMMIES ARE CONSIDERED THE SAME AS HYDROGEN.
C
c DIMENSION AA(3,3),IROW(18),ISAVE(50),XANG(50)
c COMMON/OPTEF/X(50),XNAME(50),OLDF(50),D(50),HESS(50,50),
c $ IC(50),VMODE(50),MODE,ES,NSTEP,MAXSTP,IHESS,IS,
c $ NEGREQ,DMAX,DD,CONVF,CNVFMX,CONVX,CNVXMX,IUPDAT,
c $ MXHESS,EIGMAX,EIGMIN,PRNT,NRFLAG,DUMP,IDUM
c COMMON/ZMAT/IANZ(50),IZ(50,4),BL(50),ALPHA(50),BETA(50),
c $ LBL(50),LALPHA(50),LBETA(50),NZ,NVAR
C
c DATA TOANG/0.52917706/, HARTRE/4.359814D0/
c DATA A1,A2/5.38D0,4.0D0/, B1,B2/1.35D0,1.0D0/
c DATA IROW/2*1,8*2,8*3/
c DATA AA/-0.129D0,0.186D0,0.349D0,
c $ 0.186D0,0.574D0,0.805D0,
c $ 0.349D0,0.805D0,1.094D0/
c DATA IOZMAT/507/,LZMAT/351/
C
C
C RECOVER COMMON/ZMAT/
C
c CALL TREAD(IOZMAT,IANZ,LZMAT,1,LZMAT,1,0)
C
C ALLOW FOR CHANGES IN HESS WITH BASIS
C
c CALL ILSW(2,3,IA)
c AAA=A1
c IF(IA.NE.0) AAA=A2
c BBB=B1
c IF(IA.NE.0) BBB=B2
c IF(NZ.LT.2) GO TO 50
C
C CONVERT BOND LENGTHS TO ANGSTROMS
C
c DO 15 I=1,NVAR
c XANG(I)=X(I)
c IF(NREP(I,LBL,NZ).NE.0) XANG(I)=XANG(I)*TOANG
c 15 CONTINUE
C
C GUESS DIAGONAL 2ND DERIVATIVES IN MDYNE UNITS
C
c DO 20 I=2,NZ
c IVAR=IABS(LBL(I))
c IF(IVAR.EQ.0) GO TO 20
c IATNO=IANZ(I)
c IA=IROW(IATNO)
c IF(IATNO.LT.1.OR.IATNO.GT.18) IA=1
c JATOM=IZ(I,1)
c JATNO=IANZ(JATOM)
c IB=IROW(JATNO)
c IF(JATNO.LT.1.OR.JATNO.GT.18) IB=1
c IF(IC(IVAR).EQ.0) HESS(IVAR,IVAR) = HESS(IVAR,IVAR)
c $ + AAA/( (XANG(IVAR)-AA(IA,IB))**3 )
c 20 CONTINUE
c IF(NZ.LT.3) GO TO 50
C
c DO 30 I=3,NZ
c IVAR=IABS(LALPHA(I))
c IF(IVAR.EQ.0) GO TO 30
c IF(IC(IVAR).EQ.0) HESS(IVAR,IVAR)=HESS(IVAR,IVAR)+BBB
c 30 CONTINUE
c IF(NZ.LT.4) GO TO 50
C
c DO 40 I=4,NZ
c IVAR=IABS(LBETA(I))
c IF(IVAR.EQ.0) GO TO 40
c IF(IC(IVAR).EQ.0) HESS(IVAR,IVAR)=HESS(IVAR,IVAR)+BBB
c 40 CONTINUE
C
c 50 CONTINUE
C
C CONVERT TO ATOMIC UNITS
C
c CONSTR=TOANG**2/HARTRE
c CONBND=1.0D0/HARTRE
c CALL ICLEAR(50,ISAVE)
c NSAVE=0
c DO 60 I=2,NZ
c IVAR=IABS(LBL(I))
c IF(IVAR.EQ.0.OR.NREP(IVAR,ISAVE,NSAVE).NE.0) GO TO 60
c NSAVE=NSAVE+1
c ISAVE(NSAVE)=IVAR
c IF(IC(IVAR).EQ.0) HESS(IVAR,IVAR)=HESS(IVAR,IVAR)*CONSTR
c 60 CONTINUE
C
c CALL ICLEAR(NSAVE,ISAVE)
c NSAVE=0
c DO 70 I=3,NZ
c IVAR=IABS(LALPHA(I))
c IF(IVAR.EQ.0.OR.NREP(IVAR,ISAVE,NSAVE).NE.0) GO TO 70
c NSAVE=NSAVE+1
c ISAVE(NSAVE)=IVAR
c IF(IC(IVAR).EQ.0) HESS(IVAR,IVAR)=HESS(IVAR,IVAR)*CONBND
c 70 CONTINUE
C
c DO 80 I=4,NZ
c IVAR=IABS(LBETA(I))
c IF(IVAR.EQ.0.OR.NREP(IVAR,ISAVE,NSAVE).NE.0) GO TO 80
c NSAVE=NSAVE+1
c ISAVE(NSAVE)=IVAR
c IF(IC(IVAR).EQ.0) HESS(IVAR,IVAR)=HESS(IVAR,IVAR)*CONBND
c 80 CONTINUE
C
c RETURN
c END
SUBROUTINE ESTIME(EIGVAL,FX,DD,SKAL,CHNGE,NVAR)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 EIGVAL(1),FX(1)
REAL*8 LAMDA,LAMDA0
COMMON/IO/IN,IOUT
COMMON/LAMBDA/LAMDA,LAMDA0,IT
C
C ESTIMATE THE LIKELY ENERGY CHANGE DUE TO THE STEP D
C
CHNGE=0.0D0
DO 10 I=1,NVAR
DUMMY=LAMDA
IF(I.EQ.IT) DUMMY=LAMDA0
TEMP=FX(I)*FX(I)*(DUMMY-EIGVAL(I)/2)
TEMP=TEMP/( (DUMMY-EIGVAL(I))*(DUMMY-EIGVAL(I)) )
CHNGE=CHNGE+TEMP
10 CONTINUE
C
CHNGE=CHNGE/(1+DD*DD)
C
C SCALE THE ENERGY, IN CASE THE STEP WAS SCALED
C (A SCALED ENERGY CHANGE IS LIKELY TO BE AN OVERESTIMATE)
C
CHNGE=CHNGE*SKAL
C
RETURN
END
c SUBROUTINE EXITEF(ICHAIN)
c IMPLICIT REAL*8(A-H,O-Z)
c REAL*8 GEN(47)
c LOGICAL PRNT,NRFLAG,DUMP
c COMMON/IO/IN,IOUT
c COMMON/GRDNT/ENERGY,F(50),FRCNST(1275),NVAR,IGETFC
c COMMON/OPTEF/X(50),XNAME(50),OLDF(50),D(50),HESS(50,50),
c $ IC(50),VMODE(50),MODE,ES,NSTEP,MXSTEP,IHESS,IS,
c $ NEGREQ,DMAX,DD,CONVF,CNVFMX,CONVX,CNVXMX,IUPDAT,
c $ MXHESS,EIGMAX,EIGMIN,PRNT,NRFLAG,DUMP,IDUM
c COMMON/ZSUBST/ANAMES(50),VALUES(50),INTVEC(50),FPVEC(50)
C
c DATA IRWGEN/501/,LRWGEN/47/, IGRDNT/511/,LGRDNT/1327/
c DATA IZSUBS/570/,LZSUBS/175/, IOPTEF/575/,LOPTEF/2790/
C
C LET SYNOPSIS ROUTINES KNOW IF JOB IS TERMINATING
C (BY SETTING APPROPIATE FLAG IN GEN)
C
c IF(ICHAIN.EQ.1) THEN
c CALL FILEIO(2,-IRWGEN,LRWGEN,GEN,0)
c GEN(3)=2.0D0
c CALL FILEIO(1,-IRWGEN,LRWGEN,GEN,0)
C
C COPY HESS INTO FRCNST AT SAME TIME CONTRACTING FULL MATRIX
C INTO LOWER TRIANGLE
C (THIS IS SO SYNOPSIS HAS A COPY OF THE FINAL DERIVATIVES)
C
c CALL CNTRCT(FRCNST,HESS,NVAR)
c ENDIF
C
C UPDATE THE RW FILES AND EXIT
C
c CALL DVCOPY(VALUES,X,NVAR)
c CALL TWRITE(IGRDNT,ENERGY,LGRDNT,1,LGRDNT,1,0)
c CALL TWRITE(IZSUBS,ANAMES,LZSUBS,1,LZSUBS,1,0)
c CALL TWRITE(IOPTEF,X,LOPTEF,1,LOPTEF,1,0)
C
c WRITE(IOUT,1000)
c CALL CHAINX(ICHAIN)
c RETURN
C
c1000 FORMAT(/' ',24('EF-'))
C
c END
SUBROUTINE FORMD(U,EIGVAL,FX,NVAR)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 U(50,50),EIGVAL(1),FX(1)
REAL*8 LAMDA,LAMDA0,LAMDA1,LAMDA2
LOGICAL PRNT,NRFLAG,DUMP
COMMON/IO/IN,IOUT
COMMON/OPTEF/X(50),XNAME(50),OLDF(50),D(50),HESS(50,50),
$ IC(50),VMODE(50),MODE,ES,NSTEP,MXSTEP,IHESS,IS,
$ NEGREQ,DMAX,DD,CONVF,CNVFMX,CONVX,CNVXMX,IUPDAT,
$ MXHESS,EIGMAX,EIGMIN,PRNT,NRFLAG,DUMP,IDUM
COMMON/LAMBDA/LAMDA,LAMDA0,IT
C
DATA ZERO/0.0D0/, HALF/0.5D0/, TOLL/1.0D-8/
DATA STEP/0.05D0/, BIG/1.0D+3/, MAXIT/999/
C
C
C TS SEARCH FORMS A STEP BY P-RFO THAT TRIES TO MAXIMIZE
C ALONG THE DIRECTION OF A CHOSEN HESSIAN MODE
C AND MINIMIZE ALONG ALL OTHER MODES
C MIN SEARCH FORMS A STEP BY SIMPLE RFO THAT ATTEMPTS TO
C MINIMIZE ALONG ALL HESSIAN MODES
C
C
NUMIT=0
IT=0
IF(NEGREQ.EQ.0) GO TO 10
C
C (A) MAXIMIZATION ALONG ONE OF THE HESSIAN MODES
C
IF(MODE.NE.0) THEN
CALL OVRLAP(U,NMODE,NVAR)
C
C ON RETURN FROM OVRLAP, NMODE IS THE MODE ALONG WHICH
C THE ENERGY IS TO BE MAXIMIZED
C
IF(NMODE.NE.MODE) WRITE(IOUT,1000) MODE,NMODE
MODE=NMODE
IF(PRNT) WRITE(IOUT,1010) MODE
IT=MODE
C
C IF THE MODE BEING FOLLOWED IS NOW THE LOWEST MODE,
C THEN SWITCH OFF MODE FOLLOWING
C
IF(MODE.EQ.1) THEN
MODE=0
WRITE(IOUT,1015)
ENDIF
C
ELSE
C
IF(PRNT) WRITE(IOUT,1020)
IT=1
ENDIF
C
LAMDA0=EIGVAL(IT)+DSQRT(EIGVAL(IT)**2 +4.0D0*FX(IT)**2)
LAMDA0=HALF*LAMDA0
IF(PRNT) WRITE(IOUT,1030) LAMDA0
IF(NVAR.EQ.1) GO TO 40
C
C (B) MINIMIZATION ALONG ALL OTHER MODES
C
10 CONTINUE
JT=1+IT
IF(JT.GT.2) JT=1
C
IF(PRNT.AND.NEGREQ.EQ.1) WRITE(IOUT,1040)
IF(PRNT.AND.NEGREQ.EQ.0) WRITE(IOUT,1050)
C
C SOLVE ITERATIVELY FOR LAMDA
C INITIAL GUESS FOR LAMDA IS ZERO EXCEPT NOTE THAT
C LAMDA SHOULD BE LESS THAH EIGVAL(JT)
C
LAMDA=ZERO
IF(EIGVAL(JT).LT.ZERO) THEN
LAMDA=EIGVAL(JT)-STEP
LAMDA1=EIGVAL(JT)
LAMDA2=-BIG
ENDIF
C
20 NUMIT=NUMIT+1
TEMP=ZERO
DO 25 I=1,NVAR
IF(I.EQ.IT) GO TO 25
TEMP=TEMP+(FX(I)*FX(I))/(LAMDA-EIGVAL(I))
25 CONTINUE
IF(DUMP) WRITE(IOUT,1111) LAMDA,TEMP
C
C CHECK FOR CONVERGENCE OF LAMDA
C
IF(DABS(LAMDA-TEMP).LT.TOLL) GO TO 30
C
C CHECK FOR MAXIMUM ITERATIONS EXCEEDED
C
IF(NUMIT.GT.MAXIT) GO TO 91
C
C (A) SIMPLE ITERATIVE SCHEME
C (USED IF EIGVAL(JT) > ZERO)
C
IF(EIGVAL(JT).GT.ZERO) THEN
LAMDA=TEMP
GO TO 20
C
ELSE
C
C (B) CAUTIOUS BRACKETING SCHEME
C (USED IF EIGVAL(JT) < ZERO)
C
IF(TEMP.LT.LAMDA) LAMDA1=LAMDA
IF(TEMP.GT.LAMDA) LAMDA2=LAMDA
IF(LAMDA2.GT.-BIG) LAMDA=HALF*(LAMDA1+LAMDA2)
IF(LAMDA2.EQ.-BIG) LAMDA=LAMDA-STEP
GO TO 20