-
Notifications
You must be signed in to change notification settings - Fork 0
/
prepott.f
943 lines (938 loc) · 35.5 KB
/
prepott.f
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
c***********************************************************************
SUBROUTINE PREPOTT(LNPT,IAN1,IAN2,IMN1,IMN2,NPP,VLIM,RR,VV)
c** Driver subroutine of package to generate a potential function VV(i)
c at the NPP input distances RR(i) by reading, interpolating over and
c extrapolating beyond a set of up to NPTMX read-in points.
c Based subroutine PREPOT of program LEVEL.
c====================== Version of 8 June 2007 ========================
c**** Subroutine Input:
c----------------------
c LNPT is an integer specifying the operational mode:
c * LNPT > 0 : for a new case for which all potential-defining
c parameters are read in & a description printed
c * LNPT.le.0 : if potential points are to be generated in exactly
c the same manner as on preceding call, but at
c different distances RR(I) (no reads or writes)
c IAN1 & IAN2 are the atomic numbers and IMN1 & IMN2 the mass numbers
c of atoms #1 & 2, used (if needed) to specify isotope masses for
c calculating adiabatic and/or non-adiabatic BOB correction fx.
c NPP (integer) is the number of input distances RR(i) (in Angstroms)
c at which potential values VV(i) (in cm-1) are to be generated
c VLIM (cm-1) is the absolute energy at the potential asymptote
c RR (real array) is set of NPP distances where potential calculated
c----------------------
c**** Subroutine Output:
c----------------------
c VV (real array) is the set of function values generated (in cm-1)
c NCN is an integer power defining the asymptotically-dominant
c inverse-power long-range potential tail: CNN/R**NCN
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+ Calls GENINT (which calls PLYINTRP, SPLINT & SPLINE)
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c** Set maximum array dimension for the input function values to be
c interpolated over & extrapolated beyong
INCLUDE 'arrsizes.h'
INTEGER NTPMX,VMIN,ISTATE,IDAT,K
PARAMETER (NTPMX= 1600)
INTEGER I,J,IAN1,IAN2,IMN1,IMN2,INPTS,ILR,IR2,JWR,LNPT,LWR,
1 NCN,NLIN,NPP,NROW,NTP,NUSE
REAL*8 RFACT,EFACT,RH,RMIN,VLIM,VSHIFT,VV(NPP),RR(NPP),RM2(NPP),
1 XI(NTPMX),YI(NTPMX),RWR(20),VWR(20),VWRB(3),D1V(3),D1VB(3),
2 D2V(3),CNN,RDIST(8),VDIST(8),DVDR(8),D2VDR2(8)
c
c** Save variables needed for 'subsequent' LNPT.le.0 calls
SAVE ILR,IR2,NTP,NUSE
SAVE CNN,VSHIFT,XI,YI
c
DATA VWRB/3*0.D0/,D1VB/3*0.D0/
c
WRITE(6,600) VLIM
c** For a pointwise potential (NTP > 0), now read points & parameters
c controlling how the interpolation/extrapolation is to be done.
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c** NTP : define potential by interpolation over & extrapolation
c beyond the NTP read-in turning points using subroutine GENINT.
c** If NUSE > 0 interpolate with NUSE-point piecewise polynomials
c (usually choose NUSE even, say, = 6, 8 or 10). *** If(NUSE.LE.0)
c interpolate with cubic spline instead of local polynomials.
c** If IR2 > 0 , interpolate over YI*XI**2 ; otherwise on YI itself
c This may help if interpolation has trouble on steep repulsive wall.
c** ILR specifies how to extrapolate beyond largest input distance XI(i)
c If ILR < 0 , fit last 3 points to: VLIM - A*exp(-b*(R-R0)**2)
c If ILR = 0 , fit last 3 points to: VLIM - A*R**p *exp(-b*R)
c If ILR = 1 : fit last two points to: VLIM - A/R**B .
c** If(ILR > 1) fit last turning points to: VLIM - sum{of ILR
c inverse-power terms beginning with 1/R**NCN}. *** If CNN.ne.0 ,
c leading coefficient fixed at CNN ; otherwise get it from points too.
c* Assume read-in CNN value has units: [(cm-1)(Angstroms)**'NCN'].
c* If ILR = 2 or 3 , successive higher power terms differ by 1/R**2
c* If ILR > 3 : successive higher power terms differ by factor 1/R
c-------------------------------------------------------------------
READ(5,*) NTP, NUSE, IR2, ILR, NCN, CNN
c-------------------------------------------------------------------
VMIN= 1.0d9
IF(NTP.LE.0) THEN
WRITE(6,601) !'comments to say PEC generated as analytic ...'
DO I=1,NPP,8
DO J=1,8
RDIST(J)= RR(J+I-1)
ENDDO
CALL VGENP(ISTATE,RDIST,VDIST,DVDR,D2VdR2,IDAT)
DO J=1,8
VV(J+I-1)= VDIST(J)
IF(VV(J+I-1).LT.VMIN) THEN !! locate potential minimum
VMIN= VV(J+I-1)
ENDIF
ENDDO
ENDDO
IF((I+J-1).LT.NPP) THEN
DO K= I+J-1,NPP
VV(K)= VLIM
ENDDO
ENDIF
c+++ Write for testing ++++++++++++++++++++++++++++++++++
cc REWIND(10)
WRITE(10,603) (RR(I),VV(I),I= 1,NPP,20)
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++
RETURN
ENDIF
c**** End of generation of non-standard PEC *****************************
IF(NTP.GT.NTPMX) THEN
WRITE(6,602) NTP,NTPMX
STOP
ENDIF
IF(NUSE.GT.0) WRITE(6,604) NUSE,NTP
IF(NUSE.LE.0) WRITE(6,606) NTP
IF(IR2.GT.0) WRITE(6,608)
IF((ILR.GT.1).AND.(DABS(CNN).GT.0.D0))WRITE(6,610)CNN,NCN
c** Read in turning points to be interpolated over
c** RFACT & EFACT are factors required to convert units of input turning
c points (XI,YI) to Angstroms & cm-1, respectively (may = 1.d0)
c** Turning points (XI,YI) must be ordered with increasing XI(I)
c** Energy VSHIFT (cm-1) is added to the input potential points to
c make their absolute energy consistent with VLIM (often VSHIFT=Te).
c-------------------------------------------------------------------
READ(5,*) RFACT, EFACT, VSHIFT
READ(5,*) (XI(I), YI(I), I= 1,NTP)
c-------------------------------------------------------------------
WRITE(6,612) VSHIFT, RFACT, EFACT
NROW= (NTP+2)/3
DO J= 1,NROW
IF(EFACT.LE.10.D0) THEN
WRITE(6,614) (XI(I),YI(I),I= J,NTP,NROW)
ELSE
WRITE(6,616) (XI(I),YI(I),I= J,NTP,NROW)
ENDIF
ENDDO
WRITE(6,624)
DO I= 1,NTP
YI(I)= YI(I)*EFACT+ VSHIFT
XI(I)= XI(I)*RFACT
ENDDO
IF(IR2.GT.0) THEN
DO I= 1,NTP
YI(I)= YI(I)*XI(I)**2
ENDDO
ENDIF
IF((DABS(YI(NTP)-YI(NTP-1)).LE.0).AND.
1 (XI(NTP).LT.RR(NPP))) WRITE(6,618)
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
CALL GENINT(LNPT,NPP,RR,VV,NUSE,IR2,NTP,XI,YI,VLIM,ILR,NCN,CNN)
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c IF((LNPT.GT.0).AND.(LPPOT.NE.0)) THEN
c** If desired, on the first pass (i.e. if LNPT > 0) print the potential
c RH= RR(2)-RR(1)
c INPTS= IABS(LPPOT)
c IF(LPPOT.LT.0) THEN
c** Option to write resulting function compactly to channel-8.
c RMIN= RR(1)
c NLIN= NPP/INPTS+ 1
c WRITE(8,800) NLIN,VLIM
c WRITE(8,802) (RR(I),VV(I),I= 1,NPP,INPTS)
c ELSE
c** Option to print potential & its 1-st three derivatives, the latter
c calculated by differences, assuming equally spaced RR(I) values.
c WRITE(6,620)
c NLIN= NPP/(2*INPTS)+1
c RH= INPTS*RH
c DO I= 1,NLIN
c LWR= 1+ INPTS*(I-1)
c DO J= 1,2
c JWR= LWR+(J-1)*NLIN*INPTS
c IF(JWR.LE.NPP) THEN
c RWR(J)= RR(JWR)
c VWR(J)= VV(JWR)
c D1V(J)= (VWR(J)-VWRB(J))/RH
c VWRB(J)= VWR(J)
c D2V(J)= (D1V(J)-D1VB(J))/RH
c D1VB(J)= D1V(J)
c ELSE
c RWR(J)= 0.d0
c VWR(J)= 0.d0
c ENDIF
c IF(I.LE.2) THEN
c D2V(J)= 0.d0
c IF(I.EQ.1) D1V(J)= 0.d0
c ENDIF
c ENDDO
c WRITE(6,622) (RWR(J),VWR(J),D1V(J),D2V(J),J= 1,2)
c ENDDO
c ENDIF
c ENDIF
c+++ Write for testing ++++++++++++++++++++++++++++++++++
REWIND(10)
WRITE(10,603) (RR(I),VV(I),I= 1,NPP,20)
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++
IF(LNPT.GT.0) WRITE(6,624)
RETURN
600 FORMAT(' State has energy asymptote: Y(lim)=',F12.4,'[cm-1]')
601 FORMAT(/'NTP is set less than or equal to zero, so use some analyt
1ic function')
602 FORMAT(/' **** ERROR in dimensioning of arrays required'
1 ,' by GENINT; No. input points ',I5,' > NTPMX =',I4)
603 FORMAT(5x, F12.4,f14.4)
604 FORMAT(' Perform',I3,'-point piecewise polynomial interpolation ov
1er',I5,' input points' )
606 FORMAT(' Perform cubic spline interpolation over the',I5,
1 ' input points' )
608 FORMAT(' Interpolation actually performed over modified input arra
1y: Y(I) * r(I)**2')
610 FORMAT( ' Beyond read-in points extrapolate to limiting asymptotic
1 behaviour:'/20x,'Y(r) = Y(lim) - (',D16.7,')/r**',I2)
612 FORMAT(' To make input points Y(i) consistent with Y(lim), add'
1 ,' Y(shift)=',F12.4/' Scale input points: (distance)*',
2 1PD16.9,' & (energy)*',D16.9/13x,'to get required internal unit
3s [Angstroms & cm-1 for potentials]'/
4 3(' r(i) Y(i) ')/3(3X,11('--')))
614 FORMAT((3(F13.8,F12.4)))
616 FORMAT((3(F12.6,F13.8)))
618 FORMAT(/' !!! CAUTION !!! Last two mesh point YI values are equa
1l'/17x,'so extrapolation to large r will be unreliable !!!'/)
620 FORMAT(/' Function and first 2 derivatives by differences'/
1 2(' r Y(r) d1Y/dr1 d2Y/dr2')/2(2X,19('--')))
622 FORMAT(2(0PF8.3,F11.3,1PD11.3,D10.2))
c 622 FORMAT(2(0PF7.2,F12.5,1PD11.3,D10.2))
624 FORMAT(1x,38('--'))
800 FORMAT(I7,' function values with asymptotic value:',F14.6)
802 FORMAT((1X,3(F12.8,F14.6)))
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c***********************************************************************
SUBROUTINE GENINT(LNPT,NPP,XX,YY,NUSE,IR2,NTP,XI,YI,VLIM,ILR,
1 NCN,CNN)
c** GENINT produces a smooth function YY(i) at the NPP input distances
c XX(i) by performing numerical interpolation over the range of the
c NTP input function values YI(j) at the distances XI(j), and using
c analytic functions to extrapolate beyond their range to with an
c exponential at short range and a form specified by ILR, NCN & CNN
c** ILR specifies how to extrapolate beyond largest given turning pts
c If ILR < 0 , fit last 3 points to: VLIM - A*exp(-b*(R-R0)**2)
c If ILR = 0 , fit last 3 points to: VLIM - A*R**p *exp(-b*R)
c If ILR = 1 : fit last two points to: VLIM - A/R**B .
c* If(ILR.ge.2) fit last turning points to: VLIM - sum(of ILR
c inverse-power terms beginning with 1/R**NCN). *** If CNN.ne.0 ,
c leading coefficient fixed at CNN ; otherwise get it from points too.
c* Assume read-in CNN value has units: ((cm-1)(Angstroms)**'NCN').
c If ILR = 2 or 3 , successive higher power terms differ by 1/R**2
c If ILR > 3 : this factor is 1/R .
c=== Calls subroutines PLYINTRP, SPLINT & SPLINE ==================
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER I,J,IFXCN,IDER,IR2,ILR,ISR,LNPT,MBEG,MFIN,MINNER,
1 NN,NPP,NUSE,NUST,NORD,NCN,NCN2,NCN4,NTP,
2 IMX1,NMX,JR2,JMAX,MI(10),MF(10)
REAL*8 ASR,BSR,CSR,ALR,BLR,CLR,DCSR,ADCSR,PDCSR,VRAT,
1 DX1,DX2,DX3,EX1,EX2,EX3,CNN,VLIM,X1,X2,X3,Y1,Y2,Y3,
1 XX(NPP),YY(NPP),XI(NTP),YI(NTP),XJ(20),YJ(20),DUMM(20)
c
SAVE ASR,BSR,CSR,ISR,ALR,BLR,CLR,IMX1,NMX,JR2,JMAX
c
NUST= NUSE/2
IF(NUSE.LE.0) NUST= 2
IDER= 0
c** Determine if/where need to begin extrapolation beyond input data
c XX(MI(J)) is the 1-st mesh point past turning point XI(J) .
c XX(MF(J)) is the last mesh point before turning point XI(NTP+1-J)
DO 6 J = 1,NUST
MI(J)= 1
MF(J)= 0
DO I= 1,NPP
IF(XX(I).LE.XI(J)) MI(J)= I+ 1
IF(XX(I).GE.XI(NTP+1-J)) GO TO 6
MF(J)= I
ENDDO
6 CONTINUE
IF(NUST.LT.2) THEN
MFIN= MI(1)-1
ELSE
MFIN= MI(2)-1
ENDIF
IF(LNPT.GT.0) THEN
c-----------------------------------------------------------------------
c** For a new case determine analytic functions for extrapolating beyond
c the range of the input points (if necessary) on this or later calls.
c** Try to fit three innermost turning points to V(R)=A+B*DEXP(-C*R).
c** If unsatisfactory, extrapolate inward with inverse power function
IF(IR2.LE.0) THEN
DO I= 1,4
YJ(I)= YI(I)
ENDDO
ELSE
DO I= 1,4
YJ(I)= YI(I)/XI(I)**2
ENDDO
ENDIF
X1= XI(1)
X2= XI(2)
X3= XI(3)
Y1= YJ(1)
Y2= YJ(2)
Y3= YJ(3)
IF((Y1-Y2)*(Y2-Y3).LE.0.d0) THEN
c** If 3 innermost points not monotonic, use A+B/X inward extrapoln.
ISR= 0
WRITE(6,600)
ELSE
c** Use cubic through innermost points to get initial trial exponent
c from ratio of derivatives, Y''/Y'
IDER= 2
ISR= 4
CALL PLYINTRP(XI,YJ,ISR,X2,XJ,ISR,IDER)
CSR= XJ(3)/XJ(2)
DCSR= DABS(CSR*X2)
IF(DCSR.GT.1.5D+2) THEN
c** If exponential causes overflows, use inverse power inward extrapoln.
ISR= 0
WRITE(6,602) CSR
GO TO 20
ENDIF
c** Prepare parameters for inward exponential extrapolation
VRAT= (Y3- Y2)/(Y1- Y2)
DX1= X1- X2
DX3= X3- X2
EX2= 1.D0
ADCSR= 1.d99
c** Now iterate (with actual point) to get exact exponent coefficient
DO J= 1,15
PDCSR= ADCSR
EX1= DEXP( CSR*DX1)
EX3= DEXP( CSR*DX3)
DCSR= (VRAT- (EX3- EX2)/(EX1- EX2)) /
1 ((X3*EX3- X2 - (X1*EX1- X2)*(EX3-EX2)/(EX1- EX2))/(EX1- EX2))
ADCSR= ABS(DCSR)
IF((ADCSR.GT.PDCSR).AND.(ADCSR.LT.1.d-8)) GO TO 12
IF(ADCSR.LT.1.d-12) GO TO 12
CSR= CSR+ DCSR
ENDDO
WRITE(6,604) DCSR
12 BSR= (Y1-Y2)/(EX1-EX2)
ASR= Y2-BSR*EX2
BSR= BSR*DEXP(-CSR*X2)
WRITE(6,606) X2,ASR,BSR,CSR
ENDIF
20 IF(ISR.LE.0) THEN
IF((X1*X2).LE.0.d0) THEN
c** If 1'st two mesh points of opposite sign, extrapolate linearly
ISR= -1
ASR= Y2
BSR= (Y2- Y1)/(X2- X1)
CSR= X2
WRITE(6,608) X2,ASR,BSR,CSR
ELSE
c** For inward extrapolation as inverse power through 1'st two points ..
BSR= (Y1-Y2)* X1*X2/(X2- X1)
ASR= Y1-BSR/X1
CSR= X2
WRITE(6,610) X2,ASR,BSR
ENDIF
ENDIF
ENDIF
600 FORMAT(' ** CAUTION ** Exponential inward extrapolation fails'/
1 16x,'since first 3 points not monotonic, ... so ...')
602 FORMAT(' *** CAUTION ** inward extrapolation exponent coefficient
1 C=',D12.4/10x,'could cause overflows, ... so ...')
604 FORMAT(' *** CAUTION ** after 15 tries inward extrap. exponent coe
1fft change is',1PD9.1)
606 FORMAT(' Extrapolate to X .le.',F7.4,' with'/' Y=',F13.3,
1 SP,1PD15.6,' * exp(',SS,D13.6,'*X)')
608 FORMAT(' Extrapolate to X .le.',F8.4,' with'/' Y=',F13.3,
1 SP,1PD16.7,' * [X - (',SS,F8.4,')]')
610 FORMAT(' Extrapolate to X .le.',F8.4,' with Y=',F12.3,
1 SP,1PD15.6,')/X**1')
c
IF(MFIN.GT.0) THEN
c** If needed, calculate function in inner extrapolation region
IF(ISR.GT.0) THEN
c ... either as an exponential
DO I= 1,MFIN
EX1= CSR*XX(I)
IF(DABS(EX1).GT.1.D+2) EX1= 1.D+2*DSIGN(1.d0,EX1)
YY(I)= ASR+BSR*DEXP(EX1)
ENDDO
ELSEIF(ISR.EQ.0) THEN
c ... or if that fails, as an inverse power
DO I= 1,MFIN
YY(I)= ASR+BSR/XX(I)
ENDDO
ELSEIF(ISR.LT.0) THEN
c ... or if X changes sign, extrapolate inward linearly
DO I= 1,MFIN
YY(I)= ASR+ BSR*(XX(I)- CSR)
ENDDO
ENDIF
ENDIF
c** End of inward extrapolation procedure
c-----------------------------------------------------------------------
MINNER= MFIN
IF(NUST.GT.2) THEN
c** If(NUSE.gt.5) minimize spurious behaviour by interpolating with
c order less than NUSE on intervals near inner end of range
DO J= 3,NUST
NORD= 2*(J-1)
MBEG= MI(J-1)
MFIN= MI(J)-1
IF(MFIN.GE.MBEG) THEN
DO I= MBEG,MFIN
CALL PLYINTRP(XI,YI,NTP,XX(I),DUMM,NORD,IDER)
YY(I)= DUMM(1)
ENDDO
ENDIF
ENDDO
ENDIF
c** Main interpolation step begins here
c=======================================================================
MBEG= MI(NUST)
MFIN= MF(NUST)
IF(MFIN.GE.MBEG) THEN
IF(NUSE.LE.0) THEN
c** Either ... use cubic spline for main interpolation step
CALL SPLINT(LNPT,NTP,XI,YI,MBEG,MFIN,XX,YY)
ELSE
c ... or use piecewise polynomials for main interpolation step
DO I= MBEG,MFIN
CALL PLYINTRP(XI,YI,NTP,XX(I),DUMM,NUSE,IDER)
YY(I)= DUMM(1)
ENDDO
ENDIF
ENDIF
IF(MFIN.LT.NPP) THEN
IF(NUST.LE.2) THEN
c** If(NUSE.gt.5) minimize spurious behaviour by interpolating with
c order less than NUSE on intervals near outer end of range
MBEG= MF(NUST)+1
ELSE
NN= NUST-2
DO J= 1,NN
NORD= 2*(NUST-J)
MBEG= MF(NUST-J+1)+1
MFIN= MF(NUST-J)
IF(MFIN.GE.MBEG) THEN
DO I= MBEG,MFIN
CALL PLYINTRP(XI,YI,NTP,XX(I),DUMM,NORD,IDER)
YY(I)= DUMM(1)
ENDDO
END IF
ENDDO
ENDIF
ENDIF
MBEG= MFIN+1
IF((MFIN.GT.MINNER).AND.(IR2.GT.0)) THEN
c** In (IR2.gt.0) option, now remove X**2 from the interpolated function
DO I= MINNER+1,MFIN
YY(I)= YY(I)/XX(I)**2
ENDDO
ENDIF
c** Print test of smoothness at join with analytic inward extrapolation
c IF(LNPT.GT.0) THEN
c MST= MAX0(MINNER-4,1)
c MFN= MST+8
c IF(MFN.GT.NPP) MFN= NPP
c IF(MFN.GT.MFIN) MFN= MFIN
c IF(MINNER.GT.0) WRITE(6,611) X2,((XX(I),YY(I),I= J,MFN,3),
c 1 J= MST,MST+2)
c 611 FORMAT(' Verify smoothness of inner join at X=',F9.5/
c 1 (3X,3(F10.5,G15.7)))
c ENDIF
c-----------------------------------------------------------------------
c** To extrapolate potential beyond range of given turning points ...
IF(LNPT.GT.0) THEN
c** On first entry, calculate things needed for extrapolation constants
Y1= YI(NTP)
Y2= YI(NTP-1)
Y3= YI(NTP-2)
X1= XI(NTP)
X2= XI(NTP-1)
X3= XI(NTP-2)
IF(IR2.GT.0) THEN
Y1= Y1/X1**2
Y2= Y2/X2**2
Y3= Y3/X3**2
ENDIF
ENDIF
c** Check inverse-power tail power ...
IF(NCN.LE.0) NCN= 6
IF(ILR.LT.0) THEN
IF(LNPT.GT.0) THEN
C** For ILR.lt.0 use Y = VLIM - ALR * exp[-CLR*(X - BLR)**2]
EX1= DLOG((VLIM-Y1)/(VLIM-Y2))/(X1-X2)
EX2= DLOG((VLIM-Y2)/(VLIM-Y3))/(X2-X3)
BLR= (X1+X2 - (X2+X3)*EX1/EX2)/(2.d0- 2.d0*EX1/EX2)
CLR= -EX1/(X1+X2-2.d0*BLR)
ALR= (VLIM-Y1)*DEXP(CLR*(X1-BLR)**2)
WRITE(6,614) X2,VLIM,ALR,CLR,BLR
IF(CLR.LT.0.d0) THEN
c ... but replace it by an inverse power of exponent constant negative
WRITE(6,612)
ILR= 1
GO TO 50
ENDIF
ENDIF
IF(MBEG.LE.NPP) THEN
DO I= MBEG,NPP
YY(I)= VLIM- ALR*DEXP(-CLR*(XX(I) - BLR)**2)
ENDDO
ENDIF
GO TO 90
ENDIF
IF(ILR.EQ.0) THEN
c** For ILR.le.0 use Y = VLIM - ALR * X**p * exp(-CLR*X)
IF(LNPT.GT.0) THEN
EX1= DLOG((VLIM-Y1)/(VLIM-Y2))/(X1-X2)
EX2= DLOG((VLIM-Y2)/(VLIM-Y3))/(X2-X3)
DX1= DLOG(X1/X2)/(X1-X2)
DX2= DLOG(X2/X3)/(X2-X3)
BLR= (EX1-EX2)/(DX1-DX2)
CLR= BLR*DX1- EX1
ALR= (VLIM-Y1)* DEXP(CLR*X1)/X1**BLR
WRITE(6,616) X2,VLIM,ALR,BLR,CLR
IF(CLR.LT.0.d0) THEN
c ... but replace it by an inverse power of exponent constant negative
WRITE(6,612)
ILR= 1
GO TO 50
ENDIF
ENDIF
IF(MBEG.LE.NPP) THEN
DO I= MBEG,NPP
YY(I)= VLIM- ALR*XX(I)**BLR *DEXP(-CLR*XX(I))
ENDDO
ENDIF
GO TO 90
ENDIF
50 IF(ILR.EQ.1) THEN
c** For ILR=1 , use Y = VLIM + ALR/X**BLR
IF(LNPT.GT.0) THEN
BLR= DLOG((VLIM-Y2)/(VLIM-Y1))/DLOG(X1/X2)
ALR= (Y1- VLIM)*X1**BLR
NCN= NINT(BLR)
WRITE(6,618) X2,VLIM,ALR,BLR,NCN
ENDIF
IF(MBEG.LE.NPP) THEN
DO I= MBEG,NPP
YY(I)= VLIM+ ALR/XX(I)**BLR
ENDDO
ENDIF
GO TO 90
ENDIF
c** Set constants for long-range extrapolation
IFXCN= 0
IF((CNN.GT.0.d0).OR.(CNN.LT.0.d0)) IFXCN= 1
NCN2= NCN+2
IF(ILR.EQ.2) THEN
c** For ILR=2 , use Y = VLIM - CNN/X**NCN - BSR/X**(NCN+2)
c* If CNN held fixed need ILR > 2 to prevent discontinuity
IF(LNPT.GT.0) THEN
IF(IFXCN.LE.0) THEN
CNN= ((VLIM-Y1)*X1**NCN2 -
1 (VLIM-Y2)*X2**NCN2)/(X1**2-X2**2)
ENDIF
ALR= CNN
BLR= (VLIM-Y1)*X1**NCN2 - CNN*X1**2
WRITE(6,620) X2,VLIM,CNN,NCN,BLR,NCN2
ENDIF
IF(MBEG.LE.NPP) THEN
DO I= MBEG,NPP
YY(I)= VLIM-(ALR+BLR/XX(I)**2)/XX(I)**NCN
ENDDO
ENDIF
GO TO 90
ENDIF
IF(ILR.EQ.3) THEN
c** For ILR=3 , use Y = VLIM - (CN + CN2/X**2 + CN4/X**4)/X**NCN
IF(LNPT.GT.0) THEN
NCN4= NCN+4
IF(IFXCN.GT.0) THEN
ALR= CNN
BLR= (((VLIM-Y1)*X1**NCN-ALR)*X1**4-((VLIM-Y2)
1 *X2**NCN-ALR)*X2**4)/(X1**2-X2**2)
CLR= ((VLIM-Y1)*X1**NCN-ALR-BLR/X1**2)*X1**4
ELSE
EX1= X1**2
EX2= X2**2
EX3= X3**2
DX1= (VLIM-Y1)*X1**NCN4
DX2= (VLIM-Y2)*X2**NCN4
DX3= (VLIM-Y3)*X3**NCN4
BLR= (DX1-DX2)/(EX1-EX2)
ALR= (BLR-(DX2-DX3)/(EX2-EX3))/(EX1-EX3)
BLR= BLR-ALR*(EX1+EX2)
CLR= DX1-(ALR*EX1+BLR)*EX1
ENDIF
WRITE(6,622) X2,VLIM,ALR,NCN,BLR,NCN2,CLR,NCN4
ENDIF
IF(MBEG.LE.NPP) THEN
DO I= MBEG,NPP
EX2= 1.d0/XX(I)**2
YY(I)= VLIM-(ALR+EX2*(BLR+EX2*CLR))/XX(I)**NCN
ENDDO
ENDIF
GO TO 90
ENDIF
IF(ILR.GE.4) THEN
c** For ILR.ge.4, Y = VLIM-SUM(BB(K)/X**K) , (K=NCN,NMX=NCN+ILR-1)
IF(LNPT.GT.0) THEN
IF(NCN.LE.0) NCN= 1
IMX1= ILR-1
NMX= NCN+IMX1
JR2= 0
IF(IR2.GT.0) JR2= 2
IDER= 0
JMAX= ILR
IF(IFXCN.GT.0) JMAX= IMX1
WRITE(6,624) X2,ILR,NCN,VLIM
IF(IFXCN.GT.0) WRITE(6,626) NCN,CNN
ENDIF
c** Actually extrapolate with polynomial fitted to the last JMAX
c values of (VLIM - YI(I))*XI(I)**NMX , & then convert back to YY(I).
IF(MBEG.LE.NPP) THEN
J= NTP- JMAX
DO I= 1,JMAX
J= J+1
XJ(I)= XI(J)
YJ(I)= (VLIM-YI(J)/XI(J)**JR2)*XI(J)**NMX
IF(IFXCN.GT.0) YJ(I)= YJ(I)- CNN*XI(J)**IMX1
ENDDO
DO I= MBEG,NPP
CALL PLYINTRP(XJ,YJ,JMAX,XX(I),DUMM,JMAX,IDER)
YY(I)= DUMM(1)
IF(IFXCN.GT.0) YY(I)= YY(I)+ CNN*XX(I)**IMX1
YY(I)= VLIM-YY(I)/XX(I)**NMX
ENDDO
ENDIF
ENDIF
c** Finished extrapolation section.
90 CONTINUE
c** Test smoothness at outer join to analytic extrapolation function
c IF((LNPT.GT.0).AND.(MBEG.LE.NPP)) THEN
c MST= MBEG-5
c IF(MST.LT.1) MST= 1
c MFN= MST+8
c IF(MFN.GT.NPP) MFN= NPP
c WRITE(6,627) X2,((XX(I),YY(I),I= J,MFN,3),J= MST,MST+2)
c ENDIF
c 627 FORMAT(' Verify smoothness of outer join at X=',F9.5/
c 1 (3X,3(F10.5,G15.7)))
RETURN
612 FORMAT(' *** BUT *** since exponent has positive coefficient, swi
1tch form ...')
614 FORMAT(' Function for X .GE.',F8.4,' generated as'/' Y=',
1 F12.4,' - (',1PD13.6,') * exp{-',0PF10.6,' * (r -',F9.6,')**2}')
616 FORMAT(' Function for X .GE.',F8.4,' generated as'/' Y=',
1 F12.4,' - (',1PD13.6,') * r**',0PF10.6,' * exp{-(',F11.6,'*r)}')
618 FORMAT(' Extrapolate to X .GE.',F8.4,' using'/' Y=',
1 F12.4,SP,1PD15.6,'/X**(',SS,D13.6,')] , yielding NCN=',I3)
620 FORMAT(' Extrapolate to X .GE.',F8.4,' using'/' Y=',
1 F12.4,' - [',1PD13.6,'/X**',I1,SP,D14.6,'/X**',SS,I1,']')
622 FORMAT(' Extrapolate to X .GE.',F8.4,' using'/
1 ' Y=',F12.4,' - [',1PD13.6,'/X**',I1,SP,D14.6,'/X**',
2 SS,I1,SP,D14.6,'/X**',SS,I2,']')
624 FORMAT(' Function for X .GE.',F7.3,' generated by',I3,
1 '-point inverse-power interpolation'/' with leading term 1/r**
2',I1,' relative to dissociation limit YLIM=',F11.3)
626 FORMAT(' and (dimensionless) leading coefficient fixed as C',
1 I1,'=',G15.8)
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c***********************************************************************
SUBROUTINE PLYINTRP(XI,YI,NPT,RR,C,NCFT,IDER)
c* From the NPT known mesh points (XI,YI) ,given in order of increasing
c or decreasing XI(I), select the NCFT points (XJ,YJ) surrounding the
c given point RR, and by fitting an (NCFT-1)-th degree polynomial through
c them, interpolate to find the function CC(1) and its first IDER
c derivatives (CC(I+1),I=1,IDER) evaluated at RR.
c* Adapted by R.J. Le Roy from algorithm #416,Comm.A.C.M.; 27/02/1988
c=======================================================================
INTEGER I,J,K,I1,I2,IFC,IM,IDER,J1,NH,NPT,NCFT
REAL*8 RR,XX,XI(NPT),YI(NPT),C(NCFT),XJ(20),YJ(20)
c
IF((NCFT.GT.20).OR.(NCFT.GT.NPT)) GO TO 101
NH= NCFT/2
c** First locate the known mesh points (XJ,YJ) bracketing RR
I1= 1
I2= NCFT
IF(NCFT.NE.NPT) THEN
IF(XI(NPT).LE.XI(1)) THEN
DO I= 1,NPT
IM= I
IF(XI(I).LT.RR) GO TO 20
ENDDO
ELSE
DO I= 1,NPT
IM= I
IF(XI(I).GT.RR) GO TO 20
ENDDO
ENDIF
20 I1= IM-NH
IF(I1.LE.0) I1= 1
I2= I1+NCFT-1
IF(I2.GT.NPT) THEN
I2= NPT
I1= I2-NCFT+1
ENDIF
ENDIF
J= 0
DO I= I1,I2
J= J+1
XJ(J)= XI(I)-RR
YJ(J)= YI(I)
ENDDO
c** Now determine polynomial coefficients C(I).
DO I= 2,NCFT
I1= I-1
K= I1+1
DO J= 1,I1
K= K-1
YJ(K)= (YJ(K+1)-YJ(K))/(XJ(I)-XJ(K))
ENDDO
ENDDO
C(1)= YJ(1)
DO I= 2,NCFT
XX= XJ(I)
C(I)= C(I-1)
IF(I.NE.2) THEN
I1= I-1
K= I1+1
DO J= 2,I1
K= K-1
C(K)= -XX*C(K)+C(K-1)
ENDDO
ENDIF
C(1)= YJ(I)-XX*C(1)
ENDDO
c** Finally, convert polynomial coefficients to derivatives at RR.
IFC= 1
IF(IDER.GE.NCFT) IDER= NCFT-1
IF(IDER.LE.1) GO TO 99
DO I= 2,IDER
J= I+1
IFC= IFC*I
C(J)= C(J)*IFC
ENDDO
IF(J.LT.NCFT) THEN
J1= J+1
DO I= J1,NCFT
C(I)= 0.D+0
ENDDO
ENDIF
99 RETURN
101 WRITE(6,601) NCFT,NCFT,NPT
STOP
601 FORMAT(/' *** Dimensioning ERROR in PLYINTRP : either (NCFT=',
1 I2,' .GT. 20) or (NCFT=',I2,' .GT. NPT=',I3,')')
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c**********************************************************************
SUBROUTINE SPLINT(LNPT,NTP,R1,V1,MBEG,MEND,XX,YY)
c** Subroutine to generate (if LNPT.ge.0) 4*NTP coefficients CSP(J)
c of a cubic spline passing through the NTP points (R1(J),V1(J))
c and to then calculate values of the resulting function YY(I) at the
c entering abscissae values XX(I) for I=MBEG to MEND.
c** If LNPT < 0 , generate function values at the given XX(I) using
c the coefficients CSP(J) obtained and SAVEd on a preceding call.
c** Assumes both R1(J) & XX(I) are monotonic increasing.
c+++++ Calls only subroutines SPLINE and PLYINTRP ++++++++++++++++++++++
c=======================================================================
INTEGER MAXSP
PARAMETER (MAXSP=6400)
INTEGER I,IER,I1ST,IDER,JK,K,KK,LNPT,N2,N3,NIPT,NTP,MBEG,MEND
REAL*8 EPS,R2,RI,RRR,TTMP,R1(NTP),V1(NTP),CSP(MAXSP),
1 YY(MEND),XX(MEND)
SAVE CSP
c
IF(4*NTP.GT.MAXSP) THEN
WRITE(6,602) MAXSP,NTP
STOP
ENDIF
EPS= 1.D-6*(R1(2)-R1(1))
N2= 2*NTP
N3= 3*NTP
IF(LNPT.GT.0) THEN
c** On first pass for a given data set, generate spline function
c coefficients in subroutine SPLINE
c** Start by using a cubic polynomial at each end of the range to get
c the first derivative at each end for use in defining the spline.
IDER= 1
NIPT= 4
I1ST= NTP-3
CALL PLYINTRP(R1(I1ST),V1(I1ST),NIPT,R1(NTP),CSP,NIPT,IDER)
TTMP= CSP(2)
CALL PLYINTRP(R1,V1,NIPT,R1(1),CSP,NIPT,IDER)
CSP(1)= CSP(2)
CSP(2)= TTMP
c** Now call routine to actually generate spline coefficients
ccc CALL SPLINE(R1,V1,NTP,3,CSP,MAXSP,IER)
c... using Pashov 'natural spline' with zero 2'nd derivative @ end points
CALL SPLINE(R1,V1,NTP,0,CSP,MAXSP,IER)
IF(IER .NE. 0) THEN
WRITE(6,604)
STOP
ENDIF
ENDIF
IF(MEND.LT.MBEG) GO TO 99
c** Now, use spline to generate function at desired points XX(I)
DO I= MBEG,MEND
RI= XX(I)
RRR= RI-EPS
KK= 1
c** For a monotonic increasing distance array XX(I), this statement
c speeds up the search for which set of cubic coefficients to use.
IF(I.GT.MBEG) THEN
IF(XX(I).GT.XX(I-1)) KK= JK
ENDIF
DO K= KK,NTP
JK= K
IF(R1(K).GE.RRR) GO TO 64
ENDDO
64 CONTINUE
JK= JK-1
IF(JK.LT.1) JK= 1
R2= RI-R1(JK)
YY(I)= CSP(JK)+R2*(CSP(NTP+JK)+R2*(CSP(N2+JK)+R2*CSP(N3+JK)))
ENDDO
99 RETURN
602 FORMAT(' *** ERROR in SPLINT *** Array dimension MAXSP=',I4,
1 ' cannot contain spline coefficients for NTP=',I4)
604 FORMAT(' *** ERROR in generating spline coefficients in SPLINE')
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c**********************************************************************
SUBROUTINE SPLINE(X,Y,N,IOPT,C,N4,IER)
c** Subroutine for generating cubic spline coefficients
c C(J), (J=1,N4=4*N) through the N points X(I), Y(I).
c** C(I+M*N), M=0-3 are the coefficients of order 0-3 of cubic
c polynomial expanded about X(I) so as to describe the interval:
c - X(I) to X(I+1) , if X(I) in increasing order
c - X(I-1) to X(I) , if X(I) in decreasing order.
c** IOPT indicates boundary conditions used in creating the spline .
c* If (IOPT=0) second derivatives = zero at both ends of range.
c* If (IOPT=1) 1st derivative at first point X(1) fixed at C(1),
c and 2nd derivative at X(N) = zero.
c* If (IOPT=2) 1st derivative at last point X(N) fixed at C(2),
c and 2nd derivative at X(1) = zero.
c* If (IOPT=3) constrain first derivatives at end points to have
c (read in) values C(1) at X(1) & C(2) at X(N)
c** IER is the error flag. IER=0 on return if routine successful.
c-----------------------------------------------------------------------
INTEGER I,II,IER,IOH,IOL,IOPT,J,J1,J2,J3,NER,N,N4,JMP
REAL*8 A,H,R,DY2,DYA,DYB,XB,XC,YA,YB, X(N),Y(N),C(N4)
c
JMP= 1
NER= 1000
IF(N.LE.1) GO TO 250
c** Initialization
XC= X(1)
YB= Y(1)
H= 0.D0
A= 0.D0
R= 0.D0
DYB= 0.D0
NER= 2000
c
c IOL=0 - given derivative at firstpoint
c IOH=0 - given derivative at last point
c
IOL= IOPT-1
IOH= IOPT-2
IF(IOH.EQ.1) THEN
IOL= 0
IOH= 0
ENDIF
DY2= C(2)
c
c Form the system of linear equations
c and eliminate subsequentially
c
J= 1
DO I= 1,N
J2= N+I
J3= J2+N
A= H*(2.D0-A)
DYA= DYB+H*R
IF(I.GE.N) THEN
c
c set derivative dy2 at last point
c
DYB= DY2
H= 0.D0
IF(IOH.EQ.0) GOTO 200
DYB= DYA
GOTO 220
ENDIF
J= J+JMP
XB= XC
XC= X(J)
H= XC-XB
c
c II= 0 - increasing abscissae
c II= 1 - decreasing abscissae
c
II= 0
IF(H.LT.0) II= 1
IF(H.EQ.0) GO TO 250
YA= YB
YB= Y(J)
DYB= (YB-YA)/H
IF(I.LE.1) THEN
J1= II
IF(IOL.NE.0) GO TO 220
DYA= C(1)
ENDIF
200 IF(J1.NE.II) GO TO 250
A= 1.D0/(H+H+A)
220 R= A*(DYB-DYA)
C(J3)= R
A= H*A
C(J2)= A
C(I)= DYB
ENDDO
c
c back substitution of the system of linear equations
c and computation of the other coefficients
c
A= 1.D0
J1= J3+N+II-II*N
I= N
DO IOL= 1,N
XB= X(J)
H= XC-XB
XC= XB
A= A+H
YB= R
R= C(J3)-R*C(J2)
YA= R+R
C(J3)= YA+R
C(J2)= C(I)-H*(YA+YB)
C(J1)= (YB-R)/A
C(I)= Y(J)
A= 0.D0
J= J-JMP
I= I-1
J2= J2-1
J3= J3-1
J1= J3+N+II
ENDDO
IER= 0
RETURN
250 IER= NER
RETURN
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12