-
Notifications
You must be signed in to change notification settings - Fork 1
/
main_madwave3_old.f
713 lines (594 loc) · 22.9 KB
/
main_madwave3_old.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
program madwave3
*********************************************************************
** program MadWave3 **
** version 6 , February 2018 **
** parallelized for MPI and openMP **
** last modification, January 2020 **
** references: **
** J. Chem. Phys. 107 (1997) 10085 **
** J. Chem. Phys. 109 (1998) 9391 **
** J. Chem. Phys. 123 (2005) 194309 **
** J. Chem. Phys. 125 (2006) 054102 **
** J. Phys. Chem. A 113 (2009) 14488 **
** **
** For dynamical studies of A+BC reaction dynamics **
** for: a) collisions (iphoto=0) **
** b) photoinitiated (iphoto=1) **
** c) "bright" approximation (iphoto=2) **
** and: **
** state-2-state probabilities (and amplitudes): **
** a) in product jacobi coordinates (iprod=1) **
** b) in reactant Jacobi coordinates (iprod=2) **
** c) only evaluate total reaction prob (iprod=0) **
** **
** **
** including several coupled diabatic electronic states **
** **
** Using Jacobi coordinates **
** **
** 1 **
** ^ **
** R=R_2 | **
** 2<--------------------| r=R_1 **
** | **
** | **
** 0 **
** **
** Transforming either to reactant or product Jacobi **
** to solve state-to-state probabilitites **
** **
** For a general value of Jtot in a body-fixed frame **
** with the z-axis along the R=R_2 reactant Jacobi vector. **
** **
** Input
** ------
** The progagation is done in a grid of gamma, r and R **
** and a basis for Omega **
** The angular kinetic term is evaluated **
** using a DVR method for each Omega value **
** **
** units: zots **
** Chebychev propagator in K iterations **
** using only real part of the wavepacket **
** equivalent to Gray, Balint-Kurti JCP (1998) **
** **
** this program uses the files: **
** input.dat : all input data distributed among **
** different namelist: **
** /inputgridbase/
** /inputpotmass/
** /inputprod/
** /inputprocess/ **
** /inputbnd/ --> photo initiated **
** /inputcol/ --> for collisions
** /inputtime/ **
** **
** cont.data : continuation file **
** it must contain 0 0 0 1 **
** at the beginning **
** files in directories pot,dip,func **
** generated by the preparation progeam **
** main_potini **
**
** **
** the potential matrix (for several states) is provided **
** externally as in the input of the example **
** corresponding to H3 case. **
**
** Output
** ------
** S2prod.XXXXX provides energy (cm^{-1})
** total reaction probability
** total probability (equal to 1 if converged)
**
*********************************************************************
use mod_gridYpara_01y2
use mod_pot_01y2
use mod_baseYfunciones_01y2
use mod_photoini_01y2
use mod_colini_01y2
use mod_absorcion_01y2
use mod_Hphi_01y2
use mod_flux_01y2
use mod_coortrans01_02
implicit none
include "mpif.h"
character*40 filename
integer ierr
integer :: ntimes, nloop,kminelastic
integer :: it,iit0,itinter,iloop,iloop0,indt,iloop00
integer :: i,icanp,ielec,iom,iangp,ir,ir1,ir2
real*8 :: fabs,xxx,rHnorm,autocor,autotot,time,t1,t0
real*8 :: emean,emeantot
real*8 :: xnorm1,xnorm2,xnorm1tot,xnorm2tot
*********************************************************
namelist /inputtime/ntimes, nloop,kminelastic
*********************************************************
c! Initialize MPI environment and get proc's ID and number of proc in
c! the partition.
call MPI_INIT(ierr)
! call MPI_INIT_THREAD(MPI_THREAD_FUNNELED,ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, idproc, ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr)
write(filename,'("sal."i3.3)')idproc
open(6,file=filename,status='unknown')
write(6,'(40("="),//)')
write(6,'(10x,"MadWave3 version 6 ",//)')
write(6,*)' output of proc. idproc= ',idproc,' of nproc= ',nproc
write(6,'(/,40("="),//)')
! initialization of data
call input_grid
call ini_absorcion
call paralelizacion
call pot0
! determining basis
call basis
! reactants and products functions calculation
call angular_functions
call radial_functions01_read
if(iprod.gt.0.and.npun1.ne.1)then
call product_radialf_read
endif
! reading potential
call pot2
! dimensioning vectors for propagation,
! calculating kinetic energy terms, etc
call set_vectors_Hphi
! preparing initial wave packet reading electric dipole transition
if(iphoto.eq.0)then
! initial collision wave packet
call set_colini
elseif(iphoto.eq.1.or.iphoto.eq.2)then
if(iphoto.eq.1)call read_trans_dipole
call dip_bnd
do i=1,ntotproc(idproc)
rpaqproc(i)=rpaq0(i)
enddo
else
write(6,*)' iphoto =',iphoto,' out of range'
call flush(6)
call MPI_BARRIER(MPI_COMM_WORLD, IERR)
stop
endif
! initialize total flux quantities
call ini_flux
! for flux on 02 diatomics in collisions (iphoto=0) of photodissociation (iphoto>0)
if(iprod.eq.2.and.npun1.gt.1)then
call ini_transcoor
endif
* fftw3 initialization
call difs
**>> initialating propagation
write(6,'(40("="),/,10x,"Chebysev integration",/,40("="))')
write(6,*)
write(6,*)' inputtime'
write(6,*)' ---------'
open(10,file='input.dat',status='old')
read(10,nml = inputtime)
write(6,nml = inputtime)
call flush(6)
close(10)
iit0=0
indt=0
iloop0=0
open(3,file='cont.data',status='old',err=1)
read(3,*)iit0,indt,iloop0
close(3)
write(6,*)' reading initial iit0= ',iit0,indt,iloop0
& ,' in file cont.data'
call flush(6)
1 continue
if(iit0.eq.0)then
do i=1,ntotproc(idproc)
rflanz0proc(i)=rpaqproc(i)
enddo
call difs
do i=1,ntotproc(idproc)
rpaqproc(i)=(rHpaqrec(i)
& -emindlt*rflanz0proc(i))/delta2
enddo
iit0=1
else
write(name,'("cont",i2.2,".paq")')idproc
write(6,*)' Reading wvp for continuation in file= ',name
write(6,*)' at it= ',iit0
write(6,*)
open(4,file=name,status='unknown',form='unformatted')
read(4)rflanz0proc
read(4)rpaqproc
close(4)
endif
! *>> checking average energy
emean=0.d0
xnorm1=0.d0
xnorm2=0.d0
do i=1,ntotproc(idproc)
emean=emean+rflanz0proc(i)*rHpaqrec(i)
xnorm1=xnorm1+rflanz0proc(i)**2
xnorm2=xnorm2+rHpaqrec(i)**2
enddo
call MPI_REDUCE(emean,emeantot,1,MPI_REAL8,MPI_SUM
& ,0,MPI_COMM_WORLD,ierr)
call MPI_REDUCE(xnorm1,xnorm1tot,1,MPI_REAL8,MPI_SUM
& ,0,MPI_COMM_WORLD,ierr)
call MPI_REDUCE(xnorm2,xnorm2tot,1,MPI_REAL8,MPI_SUM
& ,0,MPI_COMM_WORLD,ierr)
if(idproc.eq.0)then
write(6,*)
write(6,*)' Mean energy (eV) and norm of initial wvp= '
& ,emeantot/conve1/8065.5d0
write(6,*)
call flush(6)
emean=emeantot
endif
***********************
** Main loop in time **
***********************
write(6,1111)
iloop0=iit0/ntimes+1
time=dble(iit0)
if( iwrt_wvp == 1)then
iloop00=iloop0-1
call plot_wvp(iloop00)
endif
iloop=iloop0
it=iit0
call check(time,xnorm1tot,it,iloop)
do iloop=iloop0,nloop+iloop0
t0=MPI_Wtime()
if(idproc == 0)then
if(iphoto.ge.1.or.iprod.eq.1.or.iwrt_reac_distri.eq.1)then
write(name,'("Cvj.",i4.4)')iloop
open(15,file=name,status='new',form='unformatted')
endif
if(iprod == 2)then
write(name,'("Cvjprod.",i4.4)')iloop
open(16,file=name,status='new',form='unformatted')
endif
endif
do it = iit0+1,iit0+ntimes
time=it
call difs
do i=1,ntotproc(idproc)
call indiproc(i,icanp,ielec,iom,iangp,ir,ir1,ir2)
fabs=absfr1(ir1)*absfr2(ir2)
rHnorm=(rHpaqrec(i)-emindlt*rpaqproc(i))/delta2
xxx=(2.d0*rHnorm-fabs*rflanz0proc(i))*fabs
rflanz0proc(i)=rpaqproc(i)
rpaqproc(i)=xxx
enddo
call CoefEcheby(it)
if(npun1.gt.1)call totalflux_k(it,kminelastic)
call Cvjflux_k(it,kminelastic)
if(iprod.eq.2.and.npun1.gt.1)then
call prodpaq
call prodcvj
endif
call check(time,xnorm1tot,it,iloop)
if(idproc == 0)then
if(iphoto.ge.1.or.iprod.eq.1.or.iwrt_reac_distri.eq.1)then
write(15)it,Cvj
call flush(15)
endif
if(iprod == 2)then
write(16)it,Cvjprod
call flush(16)
endif
endif
enddo ! do it
t1=MPI_Wtime()
if(idproc.eq.0)write(6,*)' time per loop= ',t1-t0
iit0=iit0+ntimes
call flush(6)
if(npun1.gt.1)call totalflux
call wvpcont(iit0,iloop,indt)
call MPI_BARRIER(MPI_COMM_WORLD, IERR)
if(idproc.eq.0)call printprob(iit0,iloop,indt)
if( iwrt_wvp == 1)then
iloop00=iloop
call plot_wvp(iloop00)
endif
if(idproc == 0)then
if(iphoto.ge.1.or.iprod.eq.1)then
close(15)
endif
if(iprod == 2)then
close(16)
endif
endif
call MPI_BARRIER(MPI_COMM_WORLD, IERR)
if(xnorm1tot.lt.1.d-8)then
write(6,*)' norm = ',xnorm1tot,' below threshold: stop '
stop
endif
enddo ! iloop,xnorm1tot
call MPI_FINALIZE(ierr)
stop
1111 format(/,11x,'k',12x,'norm',/,30('*'),/)
end program
!----------------------------------------------------------
subroutine check(time,xnorm1tot,it,iloop)
use mod_gridYpara_01y2
use mod_pot_01y2
use mod_baseYfunciones_01y2
use mod_photoini_01y2
use mod_colini_01y2
use mod_absorcion_01y2
use mod_Hphi_01y2
implicit none
include "mpif.h"
double precision :: time,autocor,autotot
real*8 :: emean,emeantot
real*8 :: xnorm1,xnorm2,xnorm1tot
integer :: it,iloop,ierr,i
!--norm
xnorm1=0.d0
do i=1,ntotproc(idproc)
xnorm1=xnorm1+rpaqproc(i)**2
enddo
call MPI_REDUCE(xnorm1,xnorm1tot,1,MPI_REAL8,MPI_SUM
& ,0,MPI_COMM_WORLD,ierr)
call MPI_BCAST(xnorm1tot,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
if(idproc.eq.0)then
write(6,*)it,xnorm1tot
endif
if(iphoto.ge.1)then
autocor=0.d0
do i=1,ntotproc(idproc)
autocor=autocor+rpaq0(i)*rpaqproc(i)
enddo
call MPI_REDUCE(autocor,autotot,1,MPI_REAL8,MPI_SUM
& ,0,MPI_COMM_WORLD,ierr)
if(idproc.eq.0)then
write(ifileauto,*)it,autotot
call flush(ifileauto)
endif
endif ! iphoto=1
return
end subroutine check
!---------------------------------------------------------------------
subroutine CoefEcheby(ikcheb)
use mod_flux_01y2
use mod_pot_01y2, only : emindlt,delta2
use mod_gridYpara_01y2, only : hbr
use mod_flux_01y2
implicit none
complex*16 :: zexpo,zfactor
real*8 :: d,E,Es,expo,deno
integer :: ie,ikcheb
! for the representation of Energy resolved quantitites from a Chevishev expansion
! zCkcheby is the factor for Energy resolved reaction probabilities and state2state prob.
! summing on the Chebyshev iterations.
d=2.d0
if(ikcheb.eq.0)d=1.d0
zfactor=dcmplx(d*hbr/delta2,0.d0)
do ie=1,nEtot
E=etotS2(ie)
Es=(E-emindlt)/delta2
expo=-dble(ikcheb)*dacos(Es)
zexpo=cdexp(dcmplx(0.d0,expo))
deno=1.d0-Es*Es
deno=dsqrt(deno)
zCkcheby(iE)=zfactor*zexpo/dcmplx(deno,0.d0)
enddo
return
end subroutine CoefEcheby
!---------------------------------------------------------------------
subroutine wvpcont(iit0,iloop,indt)
use mod_gridYpara_01y2
use mod_pot_01y2
use mod_baseYfunciones_01y2
use mod_photoini_01y2
use mod_colini_01y2
use mod_absorcion_01y2
use mod_Hphi_01y2
use mod_flux_01y2
use mod_coortrans01_02
implicit none
real*8 :: abstot
integer :: iit0,iloop,indt
* printing only if continuation files are needed
if(ncontfile.eq.1)then
**>> continuation files
write(name,'("cont",i2.2,".paq")')idproc
open(4,file=name,status='unknown',form='unformatted')
write(4)rflanz0proc
write(4)rpaqproc
close(4)
write(name,'("cont",i2.2,".S2")')idproc
open(4,file=name,status='unknown',form='unformatted')
write(4)zS2
write(4)zCR
write(4)zCdR
close(4)
if(iprod.gt.1)then
write(name,'("cont",i2.2,".S2prod")')idproc
open(4,file=name,status='unknown',form='unformatted')
write(4)zS2prod
close(4)
endif
! if(nofe.gt.0)then
! write(name,'("cont",i2.2,".zfE")')idproc
! open(4,file=name,status='unknown',form='unformatted')
! write(4)zfE
! close(4)
! endif
if(idproc.eq.0)then
open(3,file='cont.data',status='unknown')
write(3,*)iit0,indt+1,iloop
write(3,*)abstot,nflux,ncan
close(3)
endif ! idproc.eq.0
endif !cont
return
end
!---------------------------------------------------------------------
subroutine printprob(iit0,iloop,indt)
use mod_gridYpara_01y2
use mod_pot_01y2
use mod_baseYfunciones_01y2
use mod_photoini_01y2
use mod_colini_01y2
use mod_absorcion_01y2
use mod_Hphi_01y2
use mod_flux_01y2
use mod_coortrans01_02
implicit none
include "mpif.h"
integer :: ifile,iii,ielec,ivprod,iomprod,ican,iele,j,iv
integer :: ie,iom
integer :: iit0,iloop,indt
real*8 :: S2reac,Av,S2no
real*8 :: vibprod(nviniprod:nvmaxprod)
real*8 :: rotdistri(jini:jmax)
real*8 :: S2pro(nviniprod:nvmaxprod,jiniprod:jmaxprod
& ,iomminprod:iommaxprod)
complex*16 :: zzz
* writes information at each loop in time
if(npun1.gt.1)then
write(name,'("S2prod.v",i2.2,".J",i3.3,".k",i5.5)')
& nvref,Jtot,iloop
open(20,file=name,status='unknown')
else
write(name,'("S2mat.J",i3.3,".k",i5.5)')
& Jtot,iloop
open(20,file=name,status='unknown')
endif
do ie=1,netot
****> state-2-state for reactants
rotdistri(:)=0.d0
S2reac=0.d0
do ican=1,ncan
iele=nelebas(ican)
do j=j00(ican),jmax,inc
do iv=nvini,noBCstates(j,iele)
zzz=zS2(ie,iv,j,ican)/2.d0/pi
Av=dreal(zzz*dconjg(zzz))
S2(iv,j,ican)=Av*S2factor(ie,iv,j,iele)
S2reac=S2reac+S2(iv,j,ican)
rotdistri(j)=rotdistri(j)+S2(iv,j,ican)
enddo
enddo
enddo
!****> state-2-state for products
if(npun1.gt.1)then
if(iprod.eq.1)then
S2no=0.d0
vibprod(:)=0.d0
do ican=1,ncan
iele=nelebas(ican)
do j=j00(ican),jmax,inc
do iv=nvini,min0(nvmaxprod,noBCstates(j,iele))
vibprod(iv)=vibprod(iv)+S2(iv,j,ican)
enddo
enddo
enddo
elseif(iprod.gt.1)then
S2no=0.d0
vibprod(:)=0.d0
S2pro(:,:,:)=0.d0
do iom=iomminprod,iommaxprod
do j=jiniprod,jmaxprod
do iv=nviniprod,nvmaxprod
zzz=zS2prod(ie,iv,j,iom)
Av=dreal(zzz*dconjg(zzz))*0.25d0/(pi*pi)
S2pro(iv,j,iom)=Av*S2prodfac(ie,iv,j)
S2no=S2no+S2pro(iv,j,iom)
vibprod(iv)=vibprod(iv)+S2pro(iv,j,iom)
enddo
enddo
enddo
endif
***> printing total flux to products
do iv=nviniprod,nvmaxprod
if(vibprod(iv).lt.1.d-90)vibprod(iv)=0.d0
enddo
if(S2prodtot(ie).lt.1.d-90)S2prodtot(ie)=0.d0
if(S2reac.lt.1.d-90)S2reac=0.d0
if(iprod.eq.0)then
write(20,"(501(1x,e15.7))")etotS2(ie)/conve1/8065.5d0
& ,S2prodtot(ie)*photonorm
& ,(S2prodtot(ie)+S2reac)*photonorm
& ,reacfct(ie)
elseif(iprod.eq.1)then
write(20,"(501(1x,e15.7))")etotS2(ie)/conve1/8065.5d0
& ,S2reac*photonorm
& ,(S2prodtot(ie)+S2reac)*photonorm
& ,(vibprod(iv)*photonorm,iv=nvini,min0(nvmaxprod,nvmax))
else
write(20,"(501(1x,e15.7))")etotS2(ie)/conve1/8065.5d0
& ,S2prodtot(ie)*photonorm
& ,(S2prodtot(ie)+S2reac)*photonorm
& ,S2no*photonorm
& ,(vibprod(iv)*photonorm,iv=nviniprod,nvmaxprod)
endif
else ! npun1 rotdis
write(20,"(501(1x,e15.7))")etotS2(ie)/conve1/8065.5d0
& ,S2reac*photonorm
& ,(rotdistri(j)*photonorm,j=jini,jmax)
endif
enddo ! ie=1,ne
close(20)
return
end
!---------------------------------------------------------------------
subroutine plot_wvp(iloop)
use mod_flux_01y2
use mod_pot_01y2, only : emindlt,delta2
use mod_gridYpara_01y2, only : hbr
use mod_flux_01y2
implicit none
include "mpif.h"
integer :: iloop,jelec,jr2,ir1,iang,ir2,i,icanp,ielec,iangp
integer :: nnn,ierr,ir,iom
real*8 :: r2,r1
real*8 :: fun(npun1,nangu), funtot(npun1,nangu)
* writting the wvp at each iloop
* in the same angular grid as the pes
* (i.e. each nangplot points of the angular grid)
do jelec=1,nelec
if(idproc.eq.0)then
write(name,'("wvp.elec",i2.2,".i",i3.3)')jelec,iloop
open(10,file=name,status='unknown')
endif
do jr2=1,npun2
r2=rmis2+dble(jr2-1)*ah2
if(r2.lt.absr2)then
do ir1=1,npun1
do iang=1,nangu
fun(ir1,iang)=0.d0
funtot(ir1,iang)=0.d0
enddo
enddo
do i=1,ntotproc(idproc)
call indiproc(i,icanp,ielec,iom,iangp,ir,ir1,ir2)
iang=indangreal(iangp,idproc)
if(ir2.eq.jr2.and.ielec.eq.jelec)then
fun(ir1,iang)=fun(ir1,iang)+rpaqproc(i)*rpaqproc(i)
endif
enddo
nnn=npun1*nangu
call MPI_REDUCE(fun,funtot,nnn,MPI_REAL8,MPI_SUM
& ,0,MPI_COMM_WORLD,ierr)
if(idproc.eq.0)then
do ir1=1,npun1,n1plot
r1=rmis1+dble(ir1-1)*ah1
if(r1.le.absr1)then
do iang=1,nangu
if(dabs(funtot(ir1,iang)).lt.1.d-90)then
funtot(ir1,iang)=0.d0
endif
enddo
write(10,'(500(1x,e15.7))')r1,r2
& ,(funtot(ir1,iang),iang=1,nangu,nangplot)
endif
enddo
write(10,'()')
endif
endif ! r2.le.r2abs
enddo ! ir2
if(idproc.eq.0) close(10)
enddo ! jelec
return
end subroutine plot_wvp
!---------------------------------------------------------------------