-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathart1f.f
executable file
·23173 lines (22603 loc) · 742 KB
/
art1f.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
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....................art1f.f
**************************************
*
* PROGRAM ART1.0
*
* A relativistic transport (ART) model for heavy-ion collisions
*
* sp/01/04/2002
* calculates K+K- from phi decay, dimuons from phi decay
* has finite baryon density & possibilites of varying Kaon
* in-medium mass in phiproduction-annhilation channel only.
*
*
* RELEASING DATE: JAN., 1997
***************************************
*
* Bao-An Li & Che Ming Ko
* Cyclotron Institute, Texas A&M University.
* Phone: (409) 845-1411
* e-mail: [email protected] & [email protected]
* http://wwwcyc.tamu.edu/~bali
***************************************
* Speical notice on the limitation of the code:
*
* (1) ART is a hadronic transport model
*
* (2) E_beam/A <= 15 GeV
*
* (3) The mass of the colliding system is limited by the dimensions of arrays
* which can be extended purposely. Presently the dimensions are large enough
* for running Au+Au at 15 GeV/A.
*
* (4) The production and absorption of antiparticles (e.g., ki-, anti-nucleons,
* etc) are not fully included in this version of the model. They, however,
* have essentially no effect on the reaction dynamics and observables
* related to nucleons, pions and kaons (K+) at and below AGS energies.
*
* (5) Bose enhancement for mesons and Pauli blocking for fermions are
* turned off.
*
*********************************
*
* USEFUL REFERENCES ON PHYSICS AND NUMERICS OF NUCLEAR TRANSPORT MODELS:
* G.F. BERTSCH AND DAS GUPTA, PHYS. REP. 160 (1988) 189.
* B.A. LI AND W. BAUER, PHYS. REV. C44 (1991) 450.
* B.A. LI, W. BAUER AND G.F. BERTSCH, PHYS. REV. C44 (1991) 2095.
* P. DANIELEWICZ AND G.F. BERTSCH, NUCL. PHYS. A533 (1991) 712.
*
* MAIN REFERENCES ON THIS VERSION OF ART MODEL:
* B.A. LI AND C.M. KO, PHYS. REV. C52 (1995) 2037;
* NUCL. PHYS. A601 (1996) 457.
*
**********************************
**********************************
* VARIABLES IN INPUT-SECTION: *
* *
* 1) TARGET-RELATED QUANTITIES *
* MASSTA, ZTA - TARGET MASS IN AMU, TARGET CHARGE (INTEGER) *
* *
* 2) PROJECTILE-RELATED QUANTITIES *
* MASSPR, ZPR - PROJECTILE MASS IN AMU, PROJ. CHARGE(INTEGER) *
* ELAB - BEAM ENERGY IN [MEV/NUCLEON] (REAL) *
* ZEROPT - DISPLACEMENT OF THE SYSTEM IN Z-DIREC. [FM](REAL) *
* B - IMPACT PARAMETER [FM] (REAL) *
* *
* 3) PROGRAM-CONTROL PARAMETERS *
* ISEED - SEED FOR RANDOM NUMBER GENERATOR (INTEGER) *
* DT - TIME-STEP-SIZE [FM/C] (REAL) *
* NTMAX - TOTAL NUMBER OF TIMESTEPS (INTEGER) *
* ICOLL - (= 1 -> MEAN FIELD ONLY, *
* - =-1 -> CACADE ONLY, ELSE FULL ART) (INTEGER) *
* NUM - NUMBER OF TESTPARTICLES PER NUCLEON (INTEGER) *
* INSYS - (=0 -> LAB-SYSTEM, ELSE C.M. SYSTEM) (INTEGER) *
* IPOT - 1 -> SIGMA=2; 2 -> SIGMA=4/3; 3 -> SIGMA=7/6 *
* IN MEAN FIELD POTENTIAL (INTEGER) *
* MODE - (=1 -> interpolation for pauli-blocking, *
* =2 -> local lookup, other -> unblocked)(integer) *
* DX,DY,DZ - widths of cell for paulat in coor. sp. [fm](real) *
* DPX,DPY,DPZ-widths of cell for paulat in mom. sp.[GeV/c](real) *
* IAVOID - (=1 -> AVOID FIRST COLL. WITHIN SAME NUCL. *
* =0 -> ALLOW THEM) (INTEGER) *
* IMOMEN - FLAG FOR CHOICE OF INITIAL MOMENTUM DISTRIBUTION *
* (=1 -> WOODS-SAXON DENSITY AND LOCAL THOMAS-FERMI *
* =2 -> NUCLEAR MATTER DEN. AND LOCAL THOMAS-FERMI *
* =3 -> COHERENT BOOST IN Z-DIRECTION) (INTEGER) *
* 4) CONTROL-PRINTOUT OPTIONS *
* NFREQ - NUMBER OF TIMSTEPS AFTER WHICH PRINTOUT *
* IS REQUIRED OR ON-LINE ANALYSIS IS PERFORMED *
* ICFLOW =1 PERFORM ON-LINE FLOW ANALYSIS EVERY NFREQ STEPS *
* ICRHO =1 PRINT OUT THE BARYON,PION AND ENERGY MATRIX IN *
* THE REACTION PLANE EVERY NFREQ TIME-STEPS *
* 5)
* CYCBOX - ne.0 => cyclic boundary conditions;boxsize CYCBOX *
*
**********************************
* Lables of particles used in this code *
**********************************
*
* LB(I) IS USED TO LABEL PARTICLE'S CHARGE STATE
*
* LB(I) =
clin-11/07/00:
* -30 K*-
clin-8/29/00
* -13 anti-N*(+1)(1535),s_11
* -12 anti-N*0(1535),s_11
* -11 anti-N*(+1)(1440),p_11
* -10 anti-N*0(1440), p_11
* -9 anti-DELTA+2
* -8 anti-DELTA+1
* -7 anti-DELTA0
* -6 anti-DELTA-1
clin-8/29/00-end
cbali2/7/99
* -2 antineutron
* -1 antiproton
cbali2/7/99 end
* 0 eta
* 1 PROTON
* 2 NUETRON
* 3 PION-
* 4 PION0
* 5 PION+
* 6 DELTA-1
* 7 DELTA0
* 8 DELTA+1
* 9 DELTA+2
* 10 N*0(1440), p_11
* 11 N*(+1)(1440),p_11
* 12 N*0(1535),s_11
* 13 N*(+1)(1535),s_11
* 14 LAMBDA
* 15 sigma-, since we used isospin averaged xsection for
* 16 sigma0 sigma associated K+ production, sigma0 and
* 17 sigma+ sigma+ are counted as sigma-
* 21 kaon-
* 23 KAON+
* 24 kaon0
* 25 rho-
* 26 rho0
* 27 rho+
* 28 omega meson
* 29 phi
clin-11/07/00:
* 30 K*+
* sp01/03/01
* -14 LAMBDA(bar)
* -15 sigma-(bar)
* -16 sigma0(bar)
* -17 sigma+(bar)
* 31 eta-prime
* 40 cascade-
* -40 cascade-(bar)
* 41 cascade0
* -41 cascade0(bar)
* 45 Omega baryon
* -45 Omega baryon(bar)
* sp01/03/01 end
clin-5/2008:
* 42 Deuteron (same in ampt.dat)
* -42 anti-Deuteron (same in ampt.dat)
c
* ++ ------- SEE BAO-AN LI'S NOTE BOOK
**********************************
cbz11/16/98
c PROGRAM ART
SUBROUTINE ARTMN
cbz11/16/98end
**********************************
* PARAMETERS: *
* MAXPAR - MAXIMUM NUMBER OF PARTICLES PROGRAM CAN HANDLE *
* MAXP - MAXIMUM NUMBER OF CREATED MESONS PROGRAM CAN HANDLE *
* MAXR - MAXIMUM NUMBER OF EVENTS AT EACH IMPACT PARAMETER *
* MAXX - NUMBER OF MESHPOINTS IN X AND Y DIRECTION = 2 MAXX + 1 *
* MAXZ - NUMBER OF MESHPOINTS IN Z DIRECTION = 2 MAXZ + 1 *
* AMU - 1 ATOMIC MASS UNIT "GEV/C**2" *
* MX,MY,MZ - MESH SIZES IN COORDINATE SPACE [FM] FOR PAULI LATTICE *
* MPX,MPY,MPZ- MESH SIZES IN MOMENTUM SPACE [GEV/C] FOR PAULI LATTICE *
*---------------------------------------------------------------------- *
clin PARAMETER (maxpar=200000,MAXR=50,AMU= 0.9383,
PARAMETER (MAXSTR=150001,MAXR=1,AMU= 0.9383,
1 AKA=0.498,etaM=0.5475)
PARAMETER (MAXX = 20, MAXZ = 24)
PARAMETER (ISUM = 1001, IGAM = 1100)
parameter (MX=4,MY=4,MZ=8,MPX=4,MPY=4,mpz=10,mpzp=10)
clin PARAMETER (MAXP = 14000)
*----------------------------------------------------------------------*
INTEGER OUTPAR, zta,zpr
COMMON /AA/ R(3,MAXSTR)
cc SAVE /AA/
COMMON /BB/ P(3,MAXSTR)
cc SAVE /BB/
COMMON /CC/ E(MAXSTR)
cc SAVE /CC/
COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
& RHOP(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
& RHON(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ)
cc SAVE /DD/
COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
cc SAVE /EE/
COMMON /HH/ PROPER(MAXSTR)
cc SAVE /HH/
common /ff/f(-mx:mx,-my:my,-mz:mz,-mpx:mpx,-mpy:mpy,-mpz:mpzp)
cc SAVE /ff/
common /gg/ dx,dy,dz,dpx,dpy,dpz
cc SAVE /gg/
COMMON /INPUT/ NSTAR,NDIRCT,DIR
cc SAVE /INPUT/
COMMON /PP/ PRHO(-20:20,-24:24)
COMMON /QQ/ PHRHO(-MAXZ:MAXZ,-24:24)
COMMON /RR/ MASSR(0:MAXR)
cc SAVE /RR/
common /ss/ inout(20)
cc SAVE /ss/
common /zz/ zta,zpr
cc SAVE /zz/
COMMON /RUN/ NUM
cc SAVE /RUN/
clin-4/2008:
c COMMON /KKK/ TKAON(7),EKAON(7,0:200)
COMMON /KKK/ TKAON(7),EKAON(7,0:2000)
cc SAVE /KKK/
COMMON /KAON/ AK(3,50,36),SPECK(50,36,7),MF
cc SAVE /KAON/
COMMON/TABLE/ xarray(0:1000),earray(0:1000)
cc SAVE /TABLE/
common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
cc SAVE /input1/
COMMON /DDpi/ piRHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ)
cc SAVE /DDpi/
common /tt/ PEL(-maxx:maxx,-maxx:maxx,-maxz:maxz)
&,rxy(-maxx:maxx,-maxx:maxx,-maxz:maxz)
cc SAVE /tt/
clin-4/2008:
c DIMENSION TEMP(3,MAXSTR),SKAON(7),SEKAON(7,0:200)
DIMENSION TEMP(3,MAXSTR),SKAON(7),SEKAON(7,0:2000)
cbz12/2/98
COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
& IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
cc SAVE /INPUT2/
COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX
cc SAVE /INPUT3/
cbz12/2/98end
cbz11/16/98
COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
cc SAVE /ARPRNT/
c.....note in the below, since a common block in ART is called EE,
c.....the variable EE in /ARPRC/is changed to PEAR.
clin-9/29/03 changed name in order to distinguish from /prec2/
c COMMON /ARPRC/ ITYPAR(MAXSTR),
c & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
c & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
c & XMAR(MAXSTR)
cc SAVE /ARPRC/
clin-9/29/03-end
COMMON /ARERCP/PRO1(MAXSTR, MAXR)
cc SAVE /ARERCP/
COMMON /ARERC1/MULTI1(MAXR)
cc SAVE /ARERC1/
COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
& GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
& FT1(MAXSTR, MAXR),
& PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
& EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
cc SAVE /ARPRC1/
c
DIMENSION NPI(MAXR)
DIMENSION RT(3, MAXSTR, MAXR), PT(3, MAXSTR, MAXR)
& , ET(MAXSTR, MAXR), LT(MAXSTR, MAXR), PROT(MAXSTR, MAXR)
EXTERNAL IARFLV, INVFLV
cbz11/16/98end
common /lastt/itimeh,bimp
cc SAVE /lastt/
common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG
cc SAVE /snn/
COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
cc SAVE /hbt/
common/resdcy/NSAV,iksdcy
cc SAVE /resdcy/
COMMON/RNDF77/NSEED
cc SAVE /RNDF77/
COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
clin-4/2008 zet() expanded to avoid out-of-bound errors:
real zet(-45:45)
SAVE
data zet /
4 1.,0.,0.,0.,0.,
3 1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
2 -1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1 0.,0.,0.,-1.,0.,1.,0.,-1.,0.,-1.,
s 0.,-2.,-1.,0.,1.,0.,0.,0.,0.,-1.,
e 0.,
s 1.,0.,-1.,0.,1.,-1.,0.,1.,2.,0.,
1 1.,0.,1.,0.,-1.,0.,1.,0.,0.,0.,
2 -1.,0.,1.,0.,-1.,0.,1.,0.,0.,1.,
3 0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,
4 0.,0.,0.,0.,-1./
nlast=0
do 1002 i=1,MAXSTR
ftsv(i)=0.
do 1101 irun=1,maxr
ftsvt(i,irun)=0.
1101 continue
lblast(i)=999
do 1001 j=1,4
clin-4/2008 bugs pointed out by Vander Molen & Westfall:
c xlast(i,j)=0.
c plast(i,j)=0.
xlast(j,i)=0.
plast(j,i)=0.
1001 continue
1002 continue
*-------------------------------------------------------------------*
* Input information about the reaction system and contral parameters*
*-------------------------------------------------------------------*
* input section starts here *
*-------------------------------------------------------------------*
cbz12/2/98
c.....input section is moved to subroutine ARTSET
cbz12/2/98end
*-----------------------------------------------------------------------*
* input section ends here *
*-----------------------------------------------------------------------*
* read in the table for gengrating the transverse momentum
* IN THE NN-->DDP PROCESS
call tablem
* several control parameters, keep them fixed in this code.
ikaon=1
nstar=1
ndirct=0
dir=0.02
asy=0.032
ESBIN=0.04
MF=36
*----------------------------------------------------------------------*
c CALL FRONT(12,MASSTA,MASSPR,ELAB)
*----------------------------------------------------------------------*
RADTA = 1.124 * FLOAT(MASSTA)**(1./3.)
RADPR = 1.124 * FLOAT(MASSPR)**(1./3.)
ZDIST = RADTA + RADPR
c if ( cycbox.ne.0 ) zdist=0
BMAX = RADTA + RADPR
MASS = MASSTA + MASSPR
NTOTAL = NUM * MASS
*
IF (NTOTAL .GT. MAXSTR) THEN
WRITE(12,'(//10X,''**** FATAL ERROR: TOO MANY TEST PART. ****'//
& ' '')')
STOP
END IF
*
*-----------------------------------------------------------------------
* RELATIVISTIC KINEMATICS
*
* 1) LABSYSTEM
*
ETA = FLOAT(MASSTA) * AMU
PZTA = 0.0
BETATA = 0.0
GAMMTA = 1.0
*
EPR = FLOAT(MASSPR) * (AMU + 0.001 * ELAB)
PZPR = SQRT( EPR**2 - (AMU * FLOAT(MASSPR))**2 )
BETAPR = PZPR / EPR
GAMMPR = 1.0 / SQRT( 1.0 - BETAPR**2 )
*
* BETAC AND GAMMAC OF THE C.M. OBSERVED IN THE LAB. FRAME
BETAC=(PZPR+PZTA)/(EPR+ETA)
GAMMC=1.0 / SQRT(1.-BETAC**2)
*
c WRITE(12,'(/10x,''**** KINEMATICAL PARAMETERS ****''/)')
c WRITE(12,'(10x,''1) LAB-FRAME: TARGET PROJECTILE'')')
c WRITE(12,'(10x,'' ETOTAL "GEV" '',2F11.4)') ETA, EPR
c WRITE(12,'(10x,'' P "GEV/C" '',2F11.4)') PZTA, PZPR
c WRITE(12,'(10x,'' BETA '',2F11.4)') BETATA, BETAPR
c WRITE(12,'(10x,'' GAMMA '',2F11.4)') GAMMTA, GAMMPR
IF (INSYS .NE. 0) THEN
*
* 2) C.M. SYSTEM
*
S = (EPR+ETA)**2 - PZPR**2
xx1=4.*alog(float(massta))
xx2=4.*alog(float(masspr))
xx1=exp(xx1)
xx2=exp(xx2)
PSQARE = (S**2 + (xx1+ xx2) * AMU**4
& - 2.0 * S * AMU**2 * FLOAT(MASSTA**2 + MASSPR**2)
& - 2.0 * FLOAT(MASSTA**2 * MASSPR**2) * AMU**4)
& / (4.0 * S)
*
ETA = SQRT ( PSQARE + (FLOAT(MASSTA) * AMU)**2 )
PZTA = - SQRT(PSQARE)
BETATA = PZTA / ETA
GAMMTA = 1.0 / SQRT( 1.0 - BETATA**2 )
*
EPR = SQRT ( PSQARE + (FLOAT(MASSPR) * AMU)**2 )
PZPR = SQRT(PSQARE)
BETAPR = PZPR/ EPR
GAMMPR = 1.0 / SQRT( 1.0 - BETAPR**2 )
*
c WRITE(12,'(10x,''2) C.M.-FRAME: '')')
c WRITE(12,'(10x,'' ETOTAL "GEV" '',2F11.4)') ETA, EPR
c WRITE(12,'(10x,'' P "GEV/C" '',2F11.4)') PZTA, PZPR
c WRITE(12,'(10x,'' BETA '',2F11.4)') BETATA, BETAPR
c WRITE(12,'(10x,'' GAMMA '',2F11.4)') GAMMTA, GAMMPR
c WRITE(12,'(10x,''S "GEV**2" '',F11.4)') S
c WRITE(12,'(10x,''PSQARE "GEV/C"2 '',E14.3)') PSQARE
c WRITE(12,'(/10x,''*** CALCULATION DONE IN CM-FRAME ***''/)')
ELSE
c WRITE(12,'(/10x,''*** CALCULATION DONE IN LAB-FRAME ***''/)')
END IF
* MOMENTUM PER PARTICLE
PZTA = PZTA / FLOAT(MASSTA)
PZPR = PZPR / FLOAT(MASSPR)
* total initial energy in the N-N cms frame
ECMS0=ETA+EPR
*-----------------------------------------------------------------------
*
* Start loop over many runs of different impact parameters
* IF MANYB=1, RUN AT A FIXED IMPACT PARAMETER B0, OTHERWISE GENERATE
* MINIMUM BIAS EVENTS WITHIN THE IMPACT PARAMETER RANGE OF B_MIN AND B_MAX
DO 50000 IMANY=1,MANYB
*------------------------------------------------------------------------
* Initialize the impact parameter B
if (manyb. gt.1) then
111 BX=1.0-2.0*RANART(NSEED)
BY=1.0-2.0*RANART(NSEED)
B2=BX*BX+BY*BY
IF(B2.GT.1.0) GO TO 111
B=SQRT(B2)*(BM-BI)+BI
ELSE
B=B0
ENDIF
c WRITE(12,'(///10X,''RUN NUMBER:'',I6)') IMANY
c WRITE(12,'(//10X,''IMPACT PARAMETER B FOR THIS RUN:'',
c & F9.3,'' FM''/10X,49(''*'')/)') B
*
*-----------------------------------------------------------------------
* INITIALIZATION
*1 INITIALIZATION IN ISOSPIN SPACE FOR BOTH THE PROJECTILE AND TARGET
call coulin(masspr,massta,NUM)
*2 INITIALIZATION IN PHASE SPACE FOR THE TARGET
CALL INIT(1 ,MASSTA ,NUM ,RADTA,
& B/2. ,ZEROPT+ZDIST/2. ,PZTA,
& GAMMTA ,ISEED ,MASS ,IMOMEN)
*3.1 INITIALIZATION IN PHASE SPACE FOR THE PROJECTILE
CALL INIT(1+MASSTA,MASS ,NUM ,RADPR,
& -B/2. ,ZEROPT-ZDIST/2. ,PZPR,
& GAMMPR ,ISEED ,MASS ,IMOMEN)
*3.2 OUTPAR IS THE NO. OF ESCAPED PARTICLES
OUTPAR = 0
*3.3 INITIALIZATION FOR THE NO. OF PARTICLES IN EACH SAMPLE
* THIS IS NEEDED DUE TO THE FACT THAT PIONS CAN BE PRODUCED OR ABSORBED
MASSR(0)=0
DO 1003 IR =1,NUM
MASSR(IR)=MASS
1003 CONTINUE
*3.4 INITIALIZation FOR THE KAON SPECTRUM
* CALL KSPEC0(BETAC,GAMMC)
* calculate the local baryon density matrix
CALL DENS(IPOT,MASS,NUM,OUTPAR)
*
*-----------------------------------------------------------------------
* CONTROL PRINTOUT OF INITIAL CONFIGURATION
*
* WRITE(12,'(''********** INITIAL CONFIGURATION **********''/)')
*
c print out the INITIAL density matrix in the reaction plane
c do ix=-10,10
c do iz=-10,10
c write(1053,992)ix,iz,rho(ix,0,iz)/0.168
c end do
c end do
*-----------------------------------------------------------------------
* CALCULATE MOMENTA FOR T = 0.5 * DT
* (TO OBTAIN 2ND DEGREE ACCURACY!)
* "Reference: J. AICHELIN ET AL., PHYS. REV. C31, 1730 (1985)"
*
IF (ICOLL .NE. -1) THEN
DO 700 I = 1,NTOTAL
IX = NINT( R(1,I) )
IY = NINT( R(2,I) )
IZ = NINT( R(3,I) )
clin-4/2008 check bounds:
IF(IX.GE.MAXX.OR.IY.GE.MAXX.OR.IZ.GE.MAXZ
1 .OR.IX.LE.-MAXX.OR.IY.LE.-MAXX.OR.IZ.LE.-MAXZ) goto 700
CALL GRADU(IPOT,IX,IY,IZ,GRADX,GRADY,GRADZ)
P(1,I) = P(1,I) - (0.5 * DT) * GRADX
P(2,I) = P(2,I) - (0.5 * DT) * GRADY
P(3,I) = P(3,I) - (0.5 * DT) * GRADZ
700 CONTINUE
END IF
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
*4 INITIALIZATION OF TIME-LOOP VARIABLES
*4.1 COLLISION NUMBER COUNTERS
clin 51 RCNNE = 0
RCNNE = 0
RDD = 0
RPP = 0
rppk = 0
RPN = 0
rpd = 0
RKN = 0
RNNK = 0
RDDK = 0
RNDK = 0
RCNND = 0
RCNDN = 0
RCOLL = 0
RBLOC = 0
RDIRT = 0
RDECAY = 0
RRES = 0
*4.11 KAON PRODUCTION PROBABILITY COUNTER FOR PERTURBATIVE CALCULATIONS ONLY
DO 1005 KKK=1,5
SKAON(KKK) = 0
DO 1004 IS=1,2000
SEKAON(KKK,IS)=0
1004 CONTINUE
1005 CONTINUE
*4.12 anti-proton and anti-kaon counters
pr0=0.
pr1=0.
ska0=0.
ska1=0.
* ============== LOOP OVER ALL TIME STEPS ================ *
* STARTS HERE *
* ======================================================== *
cbz11/16/98
IF (IAPAR2(1) .NE. 1) THEN
DO 1016 I = 1, MAXSTR
DO 1015 J = 1, 3
R(J, I) = 0.
P(J, I) = 0.
1015 CONTINUE
E(I) = 0.
LB(I) = 0
cbz3/25/00
ID(I)=0
c sp 12/19/00
PROPER(I) = 1.
1016 CONTINUE
MASS = 0
cbz12/22/98
c MASSR(1) = 0
c NP = 0
c NPI = 1
NP = 0
DO 1017 J = 1, NUM
MASSR(J) = 0
NPI(J) = 1
1017 CONTINUE
DO 1019 I = 1, MAXR
DO 1018 J = 1, MAXSTR
RT(1, J, I) = 0.
RT(2, J, I) = 0.
RT(3, J, I) = 0.
PT(1, J, I) = 0.
PT(2, J, I) = 0.
PT(3, J, I) = 0.
ET(J, I) = 0.
LT(J, I) = 0
c sp 12/19/00
PROT(J, I) = 1.
1018 CONTINUE
1019 CONTINUE
cbz12/22/98end
END IF
cbz11/16/98end
DO 10000 NT = 1,NTMAX
*TEMPORARY PARTICLE COUNTERS
*4.2 PION COUNTERS : LP1,LP2 AND LP3 ARE THE NO. OF P+,P0 AND P-
LP1=0
LP2=0
LP3=0
*4.3 DELTA COUNTERS : LD1,LD2,LD3 AND LD4 ARE THE NO. OF D++,D+,D0 AND D-
LD1=0
LD2=0
LD3=0
LD4=0
*4.4 N*(1440) COUNTERS : LN1 AND LN2 ARE THE NO. OF N*+ AND N*0
LN1=0
LN2=0
*4.5 N*(1535) counters
LN5=0
*4.6 ETA COUNTERS
LE=0
*4.7 KAON COUNTERS
LKAON=0
clin-11/09/00:
* KAON* COUNTERS
LKAONS=0
*-----------------------------------------------------------------------
IF (ICOLL .NE. 1) THEN
* STUDYING BINARY COLLISIONS AMONG PARTICLES DURING THIS TIME INTERVAL *
clin-10/25/02 get rid of argument usage mismatch in relcol(.nt.):
numnt=nt
CALL RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk,
& LPN,lpd,LRHO,LOMEGA,LKN,LNNK,LDDK,LNDK,LCNND,
& LCNDN,LDIRT,LDECAY,LRES,LDOU,LDDRHO,LNNRHO,
& LNNOM,numnt,ntmax,sp,akaon,sk)
c & LNNOM,NT,ntmax,sp,akaon,sk)
clin-10/25/02-end
*-----------------------------------------------------------------------
c dilepton production from Dalitz decay
c of pi0 at final time
* if(nt .eq. ntmax) call dalitz_pi(nt,ntmax)
* *
**********************************
* Lables of collision channels *
**********************************
* LCOLL - NUMBER OF COLLISIONS (INTEGER,OUTPUT) *
* LBLOC - NUMBER OF PULI-BLOCKED COLLISIONS (INTEGER,OUTPUT) *
* LCNNE - NUMBER OF ELASTIC COLLISION (INTEGER,OUTPUT) *
* LCNND - NUMBER OF N+N->N+DELTA REACTION (INTEGER,OUTPUT) *
* LCNDN - NUMBER OF N+DELTA->N+N REACTION (INTEGER,OUTPUT) *
* LDD - NUMBER OF RESONANCE+RESONANCE COLLISIONS
* LPP - NUMBER OF PION+PION elastic COLIISIONS
* lppk - number of pion(RHO,OMEGA)+pion(RHO,OMEGA)
* -->K+K- collisions
* LPN - NUMBER OF PION+N-->KAON+X
* lpd - number of pion+n-->delta+pion
* lrho - number of pion+n-->Delta+rho
* lomega - number of pion+n-->Delta+omega
* LKN - NUMBER OF KAON RESCATTERINGS
* LNNK - NUMBER OF bb-->kAON PROCESS
* LDDK - NUMBER OF DD-->KAON PROCESS
* LNDK - NUMBER OF ND-->KAON PROCESS
***********************************
* TIME-INTEGRATED COLLISIONS NUMBERS OF VARIOUS PROCESSES
RCOLL = RCOLL + FLOAT(LCOLL)/num
RBLOC = RBLOC + FLOAT(LBLOC)/num
RCNNE = RCNNE + FLOAT(LCNNE)/num
RDD = RDD + FLOAT(LDD)/num
RPP = RPP + FLOAT(LPP)/NUM
rppk =rppk + float(lppk)/num
RPN = RPN + FLOAT(LPN)/NUM
rpd =rpd + float(lpd)/num
RKN = RKN + FLOAT(LKN)/NUM
RNNK =RNNK + FLOAT(LNNK)/NUM
RDDK =RDDK + FLOAT(LDDK)/NUM
RNDK =RNDK + FLOAT(LNDK)/NUM
RCNND = RCNND + FLOAT(LCNND)/num
RCNDN = RCNDN + FLOAT(LCNDN)/num
RDIRT = RDIRT + FLOAT(LDIRT)/num
RDECAY= RDECAY+ FLOAT(LDECAY)/num
RRES = RRES + FLOAT(LRES)/num
* AVERAGE RATES OF VARIOUS COLLISIONS IN THE CURRENT TIME STEP
ADIRT=LDIRT/DT/num
ACOLL=(LCOLL-LBLOC)/DT/num
ACNND=LCNND/DT/num
ACNDN=LCNDN/DT/num
ADECAY=LDECAY/DT/num
ARES=LRES/DT/num
ADOU=LDOU/DT/NUM
ADDRHO=LDDRHO/DT/NUM
ANNRHO=LNNRHO/DT/NUM
ANNOM=LNNOM/DT/NUM
ADD=LDD/DT/num
APP=LPP/DT/num
appk=lppk/dt/num
APN=LPN/DT/num
apd=lpd/dt/num
arh=lrho/dt/num
aom=lomega/dt/num
AKN=LKN/DT/num
ANNK=LNNK/DT/num
ADDK=LDDK/DT/num
ANDK=LNDK/DT/num
* PRINT OUT THE VARIOUS COLLISION RATES
* (1)N-N COLLISIONS
c WRITE(1010,9991)NT*DT,ACNND,ADOU,ADIRT,ADDRHO,ANNRHO+ANNOM
c9991 FORMAT(6(E10.3,2X))
* (2)PION-N COLLISIONS
c WRITE(1011,'(5(E10.3,2X))')NT*DT,apd,ARH,AOM,APN
* (3)KAON PRODUCTION CHANNELS
c WRITE(1012,9993)NT*DT,ANNK,ADDK,ANDK,APN,Appk
* (4)D(N*)+D(N*) COLLISION
c WRITE(1013,'(4(E10.3,2X))')NT*DT,ADDK,ADD,ADD+ADDK
* (5)MESON+MESON
c WRITE(1014,'(4(E10.3,2X))')NT*DT,APPK,APP,APP+APPK
* (6)DECAY AND RESONANCE
c WRITE(1016,'(3(E10.3,2X))')NT*DT,ARES,ADECAY
* (7)N+D(N*)
c WRITE(1017,'(4(E10.3,2X))')NT*DT,ACNDN,ANDK,ACNDN+ANDK
c9992 FORMAT(5(E10.3,2X))
c9993 FORMAT(6(E10.3,2X))
* PRINT OUT TIME-INTEGRATED COLLISION INFORMATION
cbz12/28/98
c write(1018,'(5(e10.3,2x),/, 4(e10.3,2x))')
c & RCNNE,RCNND,RCNDN,RDIRT,rpd,
c & RDECAY,RRES,RDD,RPP
c write(1018,'(6(e10.3,2x),/, 5(e10.3,2x))')
c & NT*DT,RCNNE,RCNND,RCNDN,RDIRT,rpd,
c & NT*DT,RDECAY,RRES,RDD,RPP
cbz12/18/98end
* PRINT OUT TIME-INTEGRATED KAON MULTIPLICITIES FROM DIFFERENT CHANNELS
c WRITE(1019,'(7(E10.3,2X))')NT*DT,RNNK,RDDK,RNDK,RPN,Rppk,
c & RNNK+RDDK+RNDK+RPN+Rppk
* *
END IF
*
* UPDATE BARYON DENSITY
*
CALL DENS(IPOT,MASS,NUM,OUTPAR)
*
* UPDATE POSITIONS FOR ALL THE PARTICLES PRESENT AT THIS TIME
*
sumene=0
ISO=0
DO 201 MRUN=1,NUM
ISO=ISO+MASSR(MRUN-1)
DO 201 I0=1,MASSR(MRUN)
I =I0+ISO
ETOTAL = SQRT( E(I)**2 + P(1,I)**2 + P(2,I)**2 +P(3,I)**2 )
sumene=sumene+etotal
C for kaons, if there is a potential
C CALCULATE THE ENERGY OF THE KAON ACCORDING TO THE IMPULSE APPROXIMATION
C REFERENCE: B.A. LI AND C.M. KO, PHYS. REV. C 54 (1996) 3283.
if(kpoten.ne.0.and.lb(i).eq.23)then
den=0.
IX = NINT( R(1,I) )
IY = NINT( R(2,I) )
IZ = NINT( R(3,I) )
clin-4/2008:
c IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
c & ABS(IZ) .LT. MAXZ) den=rho(ix,iy,iz)
IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ)
2 den=rho(ix,iy,iz)
c ecor=0.1973**2*0.255*kmul*4*3.14159*(1.+0.4396/0.938)
c etotal=sqrt(etotal**2+ecor*den)
c** G.Q Li potential form with n_s = n_b and pot(n_0)=29 MeV, m^*=m
c GeV^2 fm^3
akg = 0.1727
c GeV fm^3
bkg = 0.333
rnsg = den
ecor = - akg*rnsg + (bkg*den)**2
etotal = sqrt(etotal**2 + ecor)
endif
c
if(kpoten.ne.0.and.lb(i).eq.21)then
den=0.
IX = NINT( R(1,I) )
IY = NINT( R(2,I) )
IZ = NINT( R(3,I) )
clin-4/2008:
c IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
c & ABS(IZ) .LT. MAXZ) den=rho(ix,iy,iz)
IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ)
2 den=rho(ix,iy,iz)
c* for song potential no effect on position
c** G.Q Li potential form with n_s = n_b and pot(n_0)=29 MeV, m^*=m
c GeV^2 fm^3
akg = 0.1727
c GeV fm^3
bkg = 0.333
rnsg = den
ecor = - akg*rnsg + (bkg*den)**2
etotal = sqrt(etotal**2 + ecor)
endif
c
C UPDATE POSITIONS
R(1,I) = R(1,I) + DT*P(1,I)/ETOTAL
R(2,I) = R(2,I) + DT*P(2,I)/ETOTAL
R(3,I) = R(3,I) + DT*P(3,I)/ETOTAL
c use cyclic boundary conitions
if ( cycbox.ne.0 ) then
if ( r(1,i).gt. cycbox/2 ) r(1,i)=r(1,i)-cycbox
if ( r(1,i).le.-cycbox/2 ) r(1,i)=r(1,i)+cycbox
if ( r(2,i).gt. cycbox/2 ) r(2,i)=r(2,i)-cycbox
if ( r(2,i).le.-cycbox/2 ) r(2,i)=r(2,i)+cycbox
if ( r(3,i).gt. cycbox/2 ) r(3,i)=r(3,i)-cycbox
if ( r(3,i).le.-cycbox/2 ) r(3,i)=r(3,i)+cycbox
end if
* UPDATE THE DELTA, N* AND PION COUNTERS
LB1=LB(I)
* 1. FOR DELTA++
IF(LB1.EQ.9)LD1=LD1+1
* 2. FOR DELTA+
IF(LB1.EQ.8)LD2=LD2+1
* 3. FOR DELTA0
IF(LB1.EQ.7)LD3=LD3+1
* 4. FOR DELTA-
IF(LB1.EQ.6)LD4=LD4+1
* 5. FOR N*+(1440)
IF(LB1.EQ.11)LN1=LN1+1
* 6. FOR N*0(1440)
IF(LB1.EQ.10)LN2=LN2+1
* 6.1 FOR N*(1535)
IF((LB1.EQ.13).OR.(LB1.EQ.12))LN5=LN5+1
* 6.2 FOR ETA
IF(LB1.EQ.0)LE=LE+1
* 6.3 FOR KAONS
IF(LB1.EQ.23)LKAON=LKAON+1
clin-11/09/00: FOR KAON*
IF(LB1.EQ.30)LKAONS=LKAONS+1
* UPDATE PION COUNTER
* 7. FOR PION+
IF(LB1.EQ.5)LP1=LP1+1
* 8. FOR PION0
IF(LB1.EQ.4)LP2=LP2+1
* 9. FOR PION-
IF(LB1.EQ.3)LP3=LP3+1
201 CONTINUE
LP=LP1+LP2+LP3
LD=LD1+LD2+LD3+LD4
LN=LN1+LN2
ALP=FLOAT(LP)/FLOAT(NUM)
ALD=FLOAT(LD)/FLOAT(NUM)
ALN=FLOAT(LN)/FLOAT(NUM)
ALN5=FLOAT(LN5)/FLOAT(NUM)
ATOTAL=ALP+ALD+ALN+0.5*ALN5
ALE=FLOAT(LE)/FLOAT(NUM)
ALKAON=FLOAT(LKAON)/FLOAT(NUM)
* UPDATE MOMENTUM DUE TO COULOMB INTERACTION
if (icou .eq. 1) then
* with Coulomb interaction
iso=0
do 1026 irun = 1,num
iso=iso+massr(irun-1)
do 1021 il = 1,massr(irun)
temp(1,il) = 0.
temp(2,il) = 0.
temp(3,il) = 0.
1021 continue
do 1023 il = 1, massr(irun)
i=iso+il
if (zet(lb(i)).ne.0) then
do 1022 jl = 1,il-1
j=iso+jl
if (zet(lb(j)).ne.0) then
ddx=r(1,i)-r(1,j)
ddy=r(2,i)-r(2,j)
ddz=r(3,i)-r(3,j)
rdiff = sqrt(ddx**2+ddy**2+ddz**2)
if (rdiff .le. 1.) rdiff = 1.
grp=zet(lb(i))*zet(lb(j))/rdiff**3
ddx=ddx*grp
ddy=ddy*grp
ddz=ddz*grp
temp(1,il)=temp(1,il)+ddx
temp(2,il)=temp(2,il)+ddy
temp(3,il)=temp(3,il)+ddz
temp(1,jl)=temp(1,jl)-ddx
temp(2,jl)=temp(2,jl)-ddy
temp(3,jl)=temp(3,jl)-ddz
end if
1022 continue
end if
1023 continue
do 1025 il = 1,massr(irun)
i= iso+il
if (zet(lb(i)).ne.0) then
do 1024 idir = 1,3
p(idir,i) = p(idir,i) + temp(idir,il)
& * dt * 0.00144
1024 continue
end if
1025 continue
1026 continue
end if
* In the following, we shall:
* (1) UPDATE MOMENTA DUE TO THE MEAN FIELD FOR BARYONS AND KAONS,
* (2) calculate the thermalization, temperature in a sphere of
* radius 2.0 fm AROUND THE CM
* (3) AND CALCULATE THE NUMBER OF PARTICLES IN THE HIGH DENSITY REGION
spt=0
spz=0
ncen=0
ekin=0
NLOST = 0
MEAN=0
nquark=0
nbaryn=0
csp06/18/01
rads = 2.
zras = 0.1
denst = 0.
edenst = 0.
csp06/18/01 end
DO 6000 IRUN = 1,NUM
MEAN=MEAN+MASSR(IRUN-1)
DO 5800 J = 1,MASSR(irun)
I=J+MEAN
c
csp06/18/01
radut = sqrt(r(1,i)**2+r(2,i)**2)
if( radut .le. rads )then
if( abs(r(3,i)) .le. zras*nt*dt )then
c vols = 3.14159*radut**2*abs(r(3,i)) ! cylinder pi*r^2*l
c cylinder pi*r^2*l
vols = 3.14159*rads**2*zras
engs=sqrt(p(1,i)**2+p(2,i)**2+p(3,i)**2+e(i)**2)
gammas=1.
if(e(i).ne.0.)gammas=engs/e(i)
c rho
denst = denst + 1./gammas/vols
c energy density
edenst = edenst + engs/gammas/gammas/vols
endif
endif
csp06/18/01 end
c
drr=sqrt(r(1,i)**2+r(2,i)**2+r(3,i)**2)
if(drr.le.2.0)then
spt=spt+p(1,i)**2+p(2,i)**2
spz=spz+p(3,i)**2
ncen=ncen+1
ekin=ekin+sqrt(p(1,i)**2+p(2,i)**2+p(3,i)**2+e(i)**2)-e(i)
endif
IX = NINT( R(1,I) )
IY = NINT( R(2,I) )
IZ = NINT( R(3,I) )
C calculate the No. of particles in the high density region
clin-4/2008:
c IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
c & ABS(IZ) .LT. MAXZ) THEN
IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
if(rho(ix,iy,iz)/0.168.gt.dencut)go to 5800
if((rho(ix,iy,iz)/0.168.gt.5.).and.(e(i).gt.0.9))
& nbaryn=nbaryn+1
if(pel(ix,iy,iz).gt.2.0)nquark=nquark+1
endif
c*
c If there is a kaon potential, propogating kaons
if(kpoten.ne.0.and.lb(i).eq.23)then
den=0.
clin-4/2008:
c IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
c & ABS(IZ) .LT. MAXZ)then
IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
den=rho(ix,iy,iz)
c ecor=0.1973**2*0.255*kmul*4*3.14159*(1.+0.4396/0.938)
c etotal=sqrt(P(1,i)**2+p(2,I)**2+p(3,i)**2+e(i)**2+ecor*den)
c** for G.Q Li potential form with n_s = n_b and pot(n_0)=29 MeV
c !! GeV^2 fm^3
akg = 0.1727
c !! GeV fm^3
bkg = 0.333
rnsg = den
ecor = - akg*rnsg + (bkg*den)**2
etotal = sqrt(P(1,i)**2+p(2,I)**2+p(3,i)**2+e(i)**2 + ecor)
ecor = - akg + 2.*bkg**2*den + 2.*bkg*etotal
c** G.Q. Li potential (END)
CALL GRADUK(IX,IY,IZ,GRADXk,GRADYk,GRADZk)
P(1,I) = P(1,I) - DT * GRADXk*ecor/(2.*etotal)
P(2,I) = P(2,I) - DT * GRADYk*ecor/(2.*etotal)
P(3,I) = P(3,I) - DT * GRADZk*ecor/(2.*etotal)
endif
endif
c
if(kpoten.ne.0.and.lb(i).eq.21)then
den=0.
clin-4/2008:
c IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
c & ABS(IZ) .LT. MAXZ)then
IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
den=rho(ix,iy,iz)
CALL GRADUK(IX,IY,IZ,GRADXk,GRADYk,GRADZk)
c P(1,I) = P(1,I) - DT * GRADXk*(-0.12/0.168) !! song potential
c P(2,I) = P(2,I) - DT * GRADYk*(-0.12/0.168)
c P(3,I) = P(3,I) - DT * GRADZk*(-0.12/0.168)
c** for G.Q Li potential form with n_s = n_b and pot(n_0)=29 MeV
c !! GeV^2 fm^3
akg = 0.1727
c !! GeV fm^3
bkg = 0.333
rnsg = den
ecor = - akg*rnsg + (bkg*den)**2
etotal = sqrt(P(1,i)**2+p(2,I)**2+p(3,i)**2+e(i)**2 + ecor)
ecor = - akg + 2.*bkg**2*den - 2.*bkg*etotal
P(1,I) = P(1,I) - DT * GRADXk*ecor/(2.*etotal)
P(2,I) = P(2,I) - DT * GRADYk*ecor/(2.*etotal)
P(3,I) = P(3,I) - DT * GRADZk*ecor/(2.*etotal)
c** G.Q. Li potential (END)
endif