program lpotp ******************************************************************************* c *** Paez modification for proton scattering - revised by RHL c Further revisions by M. Sagen c Minor changes to output format by L.Berge 26 Nov 88 c Now will only use Nucleons as the projectiles. All c references to the Pion code have been removed. c *** Converted to double precision by M. Sagen (August, 1986) c *** Calculation of R-matrix and T-matrix by matrix inversion with c potential in momentum space for spin 0 and 1/2, at any energy. c *** Energies, momenta, and masses in MeV units throughout program c *** Input switches - Nifty(1) thru Nifty(20) follow c c nifty(1) = -2, (p,n) cex after p,n c = 5, protons c = 9, neutron c nifty(2)*5MeV = binding energy (-5. * nifty(2) used) for 3-body energy c nifty(3) = 0, default case--takes no action. c 1, single scatt (when n(6)used else) c 2, set t=r (nonunitarize) but with ms c 3, double scattering for 1/2*1/2 c nifty(4) = 0, Regular Triplet-Singlet Mixing c 1, Triplet Singlet Mixing turned off c nifty(5) = 0, use softh(k0) c 1, sin(theat) c 2, so(on shell),no angle tf c 5, e3b energy c 6, e3b energy with folded tmatrice c 7, e3b with AAY magic vector prescription for theta c 9, soth(ko), AAY c nifty(6) = 0, nospin flip and no calc of ds and ts (0 x 0) c 1, calculate R matrix up to and include ds c 2, calculate R matrix up to and include ts c 3, include spin flip terms (0 x 1/2) c 4, no scatt (u=0) as check of wf c 5, R matrix for only ss c 6, Ustrong = 0 c 7, Imaginary part of U = 0 c 8, spin 1/2 x 1/2 with 6 spin states c nifty(7) = 0, Set H, I, J, and K NN-partial waves to zero c = 1, Retain all NN-partial waves. c nifty(8) = 0, don't read in phases, use stored t c = 1, read in phases c nifty(9) = 0, Not used by the Proton code. c nifty(10) = 0, no Coulomb amplitude c 1, Coulomb amplitude and phase change included c 2, Coulomb amplitude, but not phase change c 3, Exact coulomb, spherical app. c 4, Exact coulomb matching, realistic shape c *** Nifty(11) may be used as temporary switch c nifty(11) = 0, Schroedinger theory c = 1, Relativistic (Dirac equation) theory c (12) = 0, for off-shell GRAZ structure c = 1, for off-shell determined via f-ratios c = 2, for PEST interaction c (13) = 0, not used for proton code c (14) = 0, not used for proton code c nifty(15) = 0, divide out proton ff from He ff c 1, dont divide out proton ff c nifty(16) = 0, modified aussian rho and rho**2 & for p* z=n c nucleus will also use t(spin)*rho(matter) c 1, Wood Saxon rho,..but Gaussian rho**2 c 2, Wood Saxon rho and rho**2 c 3, special form for He or Deuteron c 4, He3..fnsp=fnmat,fpsp=0 c 5, He3..fnsp=fnmat=fpmt, fpsp=0 c 6, use fmatter*tspin for 0*.5 (non-Gaussian nucleus) c 7, uspin=(tsp*fsp+tspn*fmt) for 0*.5 c 8, use Ray's form factor for C13 c 9, use form factor of 1 (constant) c nifty(17) = 0, KMT potential c 1, Watson potential, no A/(A-1) c nifty(18) = 0, dont turn off any 1/2*1/2 amps c 1, no E terms c 2, no A+B terms c 3, no E and no A+B terms c 4, no F term in pA amplitude c nifty(19) = 0, nifty(16) controls ff's c 1, Hadjimichael impulse approximation form factors (no MEC) c 2, full Hadjimichael (i.e. with exchange currents) c nifty(20) = 0, reduced output, NO plot file to tape10 c = 1, full output, NO plot file to tape10 c = 2, reduced output, write plot file to tape10 c = 3, full output, write plot file to tape10 c nsex =0-normal/ =1-1st time thru/ =2-2nd time thru/ =3 diff cex c *** plot up lab cross sections if nang was read in as a negative number c *** U(nk,1) = Uplus for spin flip, U(nk, 2)= Uminus c *** set = for nonflip (same code for tr,ti,trc,tic,tbr,tbi) c ---------------------------------------------- implicit real*8 (a-h, o-z) character*80 SCCSID real*8 k0 , ko , ko2 , kk real*8 Mnuc , Mnuc2 , MN , Mp , Mp2 complex*16 zf c *** N.B. need to match dimensions and # of entries in data statement c *** t's enlarged 3/84 for cex (m.s.) c *** dimension of u/f depends on ngp c %%% dimension of ur,ui (ngp+1,ngp+1,6) - LB 10 Dec 88 c %%% dimension of fr,fi (2*ngp+2,2*ngp+2) - LB 10 Dec 88 c %%% dimension of den (2*ngp+2) - LB 10 Dec 88 dimension tbr(100,12), tbi(100,12), tr(100,12),ti(100,12) dimension kk(97) dimension ur(97,97,8) , ui(97,97,8) dimension fff(100,100,220,7) c dimension zf(194,194) dimension zf(66,66) dimension den(194) , gp(96) , wt(96) , sigbrn(361) dimension theta(361) , dflip(361) , dnof(361) dimension x1(361), x2(361), x11(361), x2p(361), x3(361), $ x5(361), x4(361) dimension tva(361), tvb(361), tvc(361), tvd(361), $ tve(361), tvf(361), tma(361), tmb(361), $ tmc(361), tmd(361), tme(361), tmf(361) dimension dsig(361), chart(61,61), polar(361) dimension scoul(100) dimension nifty(20) c common /nlspfl/ ur, ui common /parallel/ fff common /foptp/ zf common /kinemt/ El , Ecm , k0, pl , s , Mp2 , Mnuc , Mnuc2 common /params/ hbarc, pi, Mp, MN, nz, na, nes, nwaves common /switch/ nifty common /sizes/ achp, acmp, wsp, achn, acmn, wsn common /ranges/ rcoul, rcut common /inputs/ Tlab , b , ymin1 , ymin2 , xnang, $ kode , lxmax , ngp , nr common /Rcomn/ rr , ri , rrb , rib common /spins/ nspina , nspinb , nspinc , nspind , nspine common /Tcoul/ Trc , Tic , Tr3c , Ti3c , Tr5c , Ti5c, xgam common /Tcomn/ tr , ti , tbr , tbi c equivalence (chart(1,1),ur(1,1,1 )) equivalence (k0,ko), (ko2,k02), c line below commented out as lu's solutin to ngp prob 1/30/92 tm c $ (sigbrn(1),gp(1)), (sigbrn(33),wt(1)), $ (dsig(1),x1(1)), (polar(1),x2(1)) c *** for ibm version remove (u,h) equivalence to avoid block data c data hbarc, pi /197.3289e+0, 3.141593e+0/ ******************************************************************************** hbarc = 197.3289e+0 pi = 3.141593e+0 icex = 1 call GetIns() c nplot = 0 if (xnang .lt. 0.) nplot = 1 xnang = abs(xnang) amass = na zcm = acmp zch = achp zws = wsp nsex = 0 nifty1 = nifty(1) c *** charge exchange, do both p-He3 and n-He3, then subtract if ( nifty1 .eq. -2 ) then nsex = 1 nifty(1) = 5 write(6,700) endif 60 continue c *** cex(proton) mod (m.s.) if (nifty(1) .eq. 5) Mp = 938.2796 if (nifty(1) .eq. 9) Mp = 939.5731 c *** Mnuc here refers to nucleus mass, MN is nucleon MN = (nz * 938.2796 + (na - nz) * 939.5731)/na Mnuc = MN * amass write(6,780) amass, Mp, MN Mnuc2 = Mnuc * Mnuc Mp2 = Mp * Mp 80 continue if (nr .le. 0) write(6,1000) c *** Now k0 for non-relativistic case is the same as for relativistic, c *** only the wave equation differs El = Tlab + Mp Ecm = SQRT( Mnuc2 + Mp2 + 2. * Mnuc * El ) k0 = SQRT( (El * El - Mp2) * Mnuc2/Ecm/Ecm ) pl = SQRT( El * El - Mp2 ) write(6,960) Tlab, El, Ecm, k0, pl write(6,980) if (nr .gt. 4) write(6,990) write(6,940) write(6,950) Tlab, k0 n1 = ngp + 1 n2 = n1 + n1 a = k0 c *** special check if b=0, use different distrb of points if (b .eq. 0.) a = 200. call GAUSS2 (ngp, kode, 0.0d0, a, gp, wt) do 90 i1 = 1 , ngp kk(i1) = gp(i1) 90 continue kk(n1) = k0 ko = k0 ko2 = ko * ko emnos = sqrt(ko2 + Mnuc2) epios = sqrt(ko2 + Mp2) if (nr.eq.0) then emnos = Mnuc epios = Mp endif s = (emnos + epios)**2 c *** Now set up the denominator for either relativistic or non-relativistic c *** case call DENOM( k0, Mnuc, emnos, epios, psfac, kk, wt, den) c *** c *** set lBorn and lmax if needed c *** lclass = ko * 1.33 * (na**(.333333))/hbarc ldmax = lxmax if (lxmax .eq. 0) ldmax = 2.5 * lclass + 2.5 lBorn = 2 * lclass if (lBorn .gt. ldmax) lBorn = ldmax xgam = 0. lBorn = ldmax c if (icex .eq. 1) lbint = lBorn c if (icex .ne. 1) lBorn = lbint icex = 2 write(6,730) ldmax, lBorn, lxmax write(6,4020) c should set lborn to lmax so born approx not made,rhl 11/88 write(6,730) ldmax, lBorn, lxmax c *** Call to PRAMPS at present necessary for He-3 in order to get proper c *** initialization though it is not clear why this should be so. IF (nifty(20) .eq. 1 .OR. nifty(20) .eq. 2) then call PRAMPS(Tlab) endif c *** Stop if doing Dirac. This is temporary. if (nifty(11) .eq. 1) STOP c *** c *** include coulomb, setup - EVEN FOR NEUTRON *** zp = 1.0 c set charge to zero for neutrons or for no coulomb case if ( nifty(1) .eq. 9 .or. nifty(10).eq.0 ) zp = 0.0 xgam = zp * nz/137.036/(ko * (emnos + epios)/emnos/epios) alph = (nz - 2)/3. alph = alph/2./(2. + 3. * alph) if (na .lt. 4) alph = 0. c *** set up sigl *** rhoc = ko * rcut/hbarc lcmax = ldmax c if changed from eq to le 9/6/90 if (lcmax .le. 1) lcmax = 2 eta = xgam call SIGCL (eta, scoul, lcmax, 100) sig0 = scoul(1) c write(6,*) ' --- eta lcmax --- ',eta,lcmax c write(6,6410)(i,scoul(i),i=1,100) c6410 format(1h ,i4,e14.6) c *** c *** do loop over l for this e and state c *** do 420 ldum = 1,ldmax ld = ldum l = ldum - 1 c*************************************************************************** c declare dimensions here to pass along to optp c*************************************************************************** c %%% changed to allow 48 integration points - LB 11 Dec 88 %%% nptmax = 97 nspdim = 6 ************** tvm 10/26/92 add aovera to call optp ****************** call optp (kk, n1, ldum, ldmax, lBorn, nptmax, nspdim, aovera) write(6,1070) l c *** c ---------- do loop over 2 - 6 spin states(if necessary)-------- c *** nspinm = 2 if (nifty(6) .eq. 8) nspinm = 8 do 410 nspin = 1,nspinm n = nspin if (n .eq. 3 .OR. n .eq. 5) go to 410 c *********** tvm 9/2/92 TS-mixing ******************************* if (n .eq. 1 .or. n .eq. 7) go to 410 nspine = nspin c *** set spin indices needed for 1/2 x 1/2 if (nifty(6) .eq. 8) then nspina = nspin - 1 nspinb = nspin nspinc = nspina + 2 * (5 - nspin) nspind = nspinc + 1 endif c *** special storage r0(11)-nspin=5 unreachable with our l sum scheme c *** store in ro(00)-nspin=2 for ldum=1(l=0) 228 if (nspin .eq. 2 $ .AND. $ nifty(6) .ne. 3 $ .AND. $ nifty(6) .ne. 8 ) go to 380 write(6,2211) n write(6,1979) ur(n1,n1,n), ui(n1,n1,n) c *** c *** test if l large enough for born approx c *** if (ldum .gt. lBorn) go to 340 c nmax = n1 c if (nspin .gt. 2) nmax = n2 c*********** TS-mixing tvm 9/8/92 c 2 above, add below nmax = n2 c *** Set up the f-matrix *** ******************* New T-Mat 8/1/92 tvm ***************** **************Add psfac to variables for fmatrx *********** call Fmatrx( den , n1 , nmax , nspin , nptmax, ldum, psfac) c *** Calculate the on-shell T matrix from inverse of f * u *** call Tmatrx ( psfac , n1 , n2 , nmax , nspin ) *^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ******************** New T-Mat 7/27/92 tvm ********************* tr(ld,nspina) = rr ti(ld,nspina) = ri tr(ld,nspinb) = rrb ti(ld,nspinb) = rib *________________________________________________________________ c *** calculation of t-matrix (on - shell), c *** t normalized to exp(i del)*sin(del) 2424 continue call coulomb( scoul, sigl, aovera, ld, nspin, njump ) 300 continue if (ldum .gt. lBorn) go to 390 c *** c ---------------born approximation ------------------ c *** 340 continue c write(6,1979) ur(n1,n1,n), ui(n1,n1,n) c call TBornA( aovera, psfac, ld, ldum, lBorn, nspin, n1 ) if (ldum .le. lBorn) go to 390 c go to 2424 c *** no spin dependence, set tl+=tl- *********** below only for 0 x 1/2 ******************* 380 continue c380 tr(ld,2) = tr(ld,1) c ti(ld,2) = ti(ld,1) c tbr(ld,2) = tbr(ld,1) c tbi(ld,2) = tbi(ld,1) c *** charge exchange check 390 if (nsex .ne. 0) then call Tcex( tr, ti, tbr, tbi, sigl, ld, nsex, n ) endif 410 continue c *** End of do loop over spin *** 420 continue c *** End of do loop over partial waves(after remove coulomb) *** c430 if (nsex .eq. 0) go to 435 c print out t matrices in every case 430 continue write(6,740) nsex do 434 ll = 1,ldmax c ************* tvm 9/2/92 TS mixing Why fight? 2 - 7 interchange in nspin *************** convention from this point on. ********************** tmdum = tr(ll,2) tvdum = ti(ll,2) tr(ll,2) = tr(ll, 7) ti(ll,2) = ti(ll, 7) tr(ll,7) = tmdum ti(ll,7) = tvdum ********* Fix iteration convention in line below, -1 tvm 7/28/92 **** 434 write(6,1031) ll-1, (tr(ll,np2), np2=1,8), (ti(ll,np2), np2=1,8) 435 continue if (nifty(6) .ne. 8) write(6,1090) if (nifty(6) .eq. 8) write(6,1091) if (nifty(6) .eq. 3) write (6,1110) c write(7, 7213) c7213 format(1h1,'th_c.m.', 4x, 'Aoolm', 7x, 'Aooml', c $ 7x, 'Aoomm', 7x, 'Aooll') c-----calculate angular disrtbn --------------------------------------- call xsects( theta, dsig, sigbrn, dnof, dflip, polar, $ scoul, alph, xgam, sig0, ldmax, noangs, $ nspin, nsex ) if (nplot .ne. 1) go to 581 c *** plot up cross sections on semilog 5 cycle paper 561 nspin = 0 if (nifty(6) .eq. 3 .OR. nifty(6) .eq. 8) nspin = 1 c *** variables for plotting titles paramet = 0.0 paramet = paramet + ((dsig(11) - 1101.0)/5.111)**2.0 * 1.0 paramet = paramet + ((dsig(13) - 897.1)/4.617)**2.0 * 1.0 paramet = paramet + ((dsig(14) - 785.4)/4.319)**2.0 paramet = paramet + ((dsig(17) - 548.4)/4.271)**2.0 * 1.0 paramet = paramet + ((dsig(23) - 183.8)/1.453)**2.0 paramet = paramet + ((dsig(25) - 115.0)/0.845)**2.0 paramet = paramet + ((dsig(35) - 1.985)/0.0237)**2.0 * 1.0 paramet = paramet + ((dsig(38) - 0.5378)/0.00961)**2.0 paramet = paramet + ((dsig(41) - 1.215)/0.0134)**2.0 * 1.0 paramet = paramet + ((dsig(51) - 1.761)/0.0128)**2.0 paramet = paramet + ((dsig(53) - 1.447)/0.0117)**2.0 paramet = paramet + ((dsig(64) - 0.103)/0.00153)**2.0 paramet = paramet + ((dsig(67) - 0.0374)/0.00078)**2.0 paramet = paramet + ((dsig(72) - 0.0289)/0.000835)**2.0 goodness = paramet / 12.0 c write(7,1139) achp, acmn, paramet, goodness c1139 format(4e12.5) c write(7,1138) paramet c1138 format(' ',e18.9) chart(1,1) = Tlab chart(2,1) = achp chart(3,1) = acmp chart(4,1) = wsp chart(5,1) = achn chart(6,1) = acmn chart(7,1) = wsn chart(8,1) = nz chart(9,1) = amass do 570 ii = 1,20 j = ii + 9 chart(j,1) = nifty(ii) 570 continue if (nifty(6) .ne. 8) then call logplt (theta, dsig, dnof, noangs, nspin, chart, ymin1) endif if (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) then call logplt (theta, dnof, dflip, noangs, 1, chart, ymin2) endif if (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) go to 581 c dsig/domega call logplt (theta, dsig, dnof, noangs, 0, chart, ymin1) call logplt (theta, dsig, dnof, noangs, 0, chart, ymin2) if(nifty(20) .lt.2) go to 576 c read in plot data from tape10 rewind 10 do 575 i = 1, noangs read(10,*) theta(i), x1(i), x2(i), x11(i), x2p(i), x3(i), $ x4(i), x5(i) 575 continue 576 continue c 2,2' 11 if(na.eq.13) go to 577 call logplt (theta, polar, dnof, noangs, -1, chart, ymin2) call logplt (theta, dflip, dnof, noangs, -1, chart, ymin2) call logplt (theta, dnof, dnof, noangs, -1, chart, ymin2) 577 continue c rhl add for c13 if(na.ne.13) go to 579 do 578 i=1,noangs 578 theta(i)= 3.*theta(i) call logplt (theta, dsig, dnof, noangs, 0, chart, ymin1) if (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) then call logplt (theta, dnof, dflip, noangs, 0, chart, ymin2) endif if (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) go to 581 call logplt (theta, dsig, dflip, noangs, 1, chart, ymin2) call logplt (theta, polar, dflip, noangs, -1, chart, ymin2) if (nifty(20).lt. 2) go to 579 c plot up plotdata stored on tape10 c call logplt (theta, x2, dflip, noangs, -1, chart, ymin2) call logplt (theta, x11, dflip, noangs, -1, chart, ymin2) call logplt (theta, x2p, dflip, noangs, -1, chart, ymin2) call logplt (theta, x3, dflip, noangs, -1, chart, ymin2) call logplt (theta, x4, dflip, noangs, -1, chart, ymin2) call logplt (theta, x5, dflip, noangs, -1, chart, ymin2) 579 continue rewind 11 do 580 i = 1, noangs c read(11,*)theta(i), tva(i), tma(i), tvb(i), tmb(i), tvc(i), c $ tmc(i), tvd(i), tmd(i), tve(i), tme(i), tvf(i), tmf(i) 580 continue c call logplt (theta, tva, tma, noangs, -2, chart, ymin2) c call logplt (theta, tvb, tmb, noangs, -2, chart, ymin2) c call logplt (theta, tvc, tmc, noangs, -2, chart, ymin2) c call logplt (theta, tvd, tmd, noangs, -2, chart, ymin2) c call logplt (theta, tve, tme, noangs, -2, chart, ymin2) c call logplt (theta, tvf, tmf, noangs, -2, chart, ymin2) 581 if (nsex .eq. 0) go to 620 if (nsex .eq. 3) go to 620 if (nsex .eq. 2) go to 610 c *** nsex = 1 pass done, setup 2nd pass c *** proton done, now do neutron nsex = 2 nifty(1) = 9 c______________________________________________________________________________ c Now reinitialize t's by calling TNCM for charge exchange c______________________________________________________________________________ energy=-200.e+00 call tncm (energy) c *** go to 60 c *** nsex = 2, two passes for cex done,now do last 610 nsex = 3 ymin1 = ymin1 - 1 go to 430 620 continue write(6,841) c stop c _______________ formats c *** c227 format(i3, 8e15.6 / 200(3x, 8e15.6 / )) 700 format (15h p, n, then cex) 730 format (1h0,'lmax(calculated) = ',i4,5x,'lBorn = ',i4,5x, $ 'lmax(read) = ',i4) 740 format (1h1,5i4) 780 format (1h0,' A=',f10.4,' Mp=',f10.4,' m(nucleon)=',f10.4//7x, $ ' Tlab', 9x, 'Eplab ', 8x, 'Ecm(sqrts)', 4x, $ 'kcm', 11x, 'plab') c820 format (47h R matrix calculated with ms series up to order,i3) 841 format(12h eof in main ,i3) c928 format(1h0,33h re/im f matrix before inversion ) c930 format (1h ,26h normalized r(k,k..k)= ,24x,2e14.7) 940 format (1h0,10x,8h Energy,14x,15h Momentum ,/) 950 format (1h ,7x,f14.4,10x,f14.4) 960 format (1h ,5f14.4) 980 format (1h0,31h relativistic calculation) 990 format (1h ,10x,33happroximate klein gordon equation) 1000 format (1h0,47h nonrelativstc calc--but with same k0 as reltv,40x 1,32h i.e. only the wave eqtndiffers) c1010 format (1h+,85x,4heta=,f10.4,8h delta=,f10.4) c1020 format (41h0********** eta greater than 1 **********,//) c1030 format (1h ,15h t(k,k,k)= ,35x,4e14.7) 1031 format (1h ,16h t(k,k,k),for l=,i3,/,10(8e10.3,/,8e10.3,/)) c1040 format (38h now with exact coulomb after matching) c1050 format (20x,40hcoulomb phase modified nuclear amplitude/16h t(k, c 1k.k) = ,4e14.7) 1070 format (1h0,30h----- angular momentum = ,i4,5h-----) c1080 format (1h ,8x,23h born approximation ) 1090 format (17h1th-c.m. cos(cm),3x,10hdsig/dw-cm,2x,11hds/dw(b.a.), 13x,1ht,4x,6hth-lab,2x,8hcos(lab),2x,5hds/dw,4x,11hds/dw(b.a.),2x, 210href(th)-cm,1x,3himf,4x,6hsignuc,4x,7hsigcoul//) c1091 format(1h1,' th-c.m. cos(cm),3x,4hq2fm,7x,10hdsig/dw-cm , c 1 1x,10hds/dw(nuc) ,1x,11hds/dw(b.a.) ,2x,5hpolar ,4x, c 2 7hnoflip ,4x,4hflip ,5x,4hpnoo ,4x,5hkolmo ,4x, c 35hclmoo ,2x,7hsigcoul /1h ,18x, c 420hza ,zb ,zc ,zd ,ze /1h ,18x, c 520hzab,zbb,zcb,zdb,zeb //) 1091 format(1h1,'th-c.m. q ',4x,'ds/dw',6x,'P=Aoono',6x, $ 'Dnono',7x,'Dloso',7x,'Dsoso',7x,'Dlolo',7x,'Aoonn',7x,'Aooon') c1093 format(1h ,f7.0,2f8.3,3e11.3,f8.3,2e11.3,3f9.3,e11.3) c1100 format(1h ,f7.0,f8.3,2e11.3,2f8.1,f10.3,2e11.3,2f9.3,e12.4 c 1,e10.2) 1110 format (1h+,21x,6hnoflip,7x,8hnoflip-b,5x,5hpolar,9x,4hflip,9x, 18hflip-b,r,//) c1120 format (1h ,18x,10e11.3) c1498 format(' ???main ur(lid,lid,1)= ',e10.4,' ui(same)= 'e10.4) c1921 format(1h ,7h sigot=,e13.5,7h sigit=,e13.5, c 27h sig2t=,e13.5,2he=,f8.2,5hxgam=,e13.5) 1979 format(1h ,17h ur ui(k0,k0,k0)= ,2e13.5) 2211 format(1h ,' .......... nspin= ',i4,' ...... ') c3335 format(1h ,19h tr ti (nspin=3-6)=,8e13.5) c4000 format(' ', f7.2, 3e11.3) c4010 format(1h ,f7.2,2e11.3) 4020 format(1h2) end