latest versions of GEMINI are available from http://wunmr.wustl.edu/~rc/ change log of fortran95 version ************************************************************** 2- May-2205 Made some minor modifications so that the progrom would compile on the Nag F95 compiler. 7- Jan-2003 bug fixed. With imf_option=2 and tl_iwbc=.true. if a fragement has small Z, then the evap.f90 routines are used rather than the light_particle subroutines - Had to modify evap.f90 so the angle calculation is sphere_lp didn't crash. 26-nov-2002 minor bug in subroutine input about polarization 2-June-2002 found bug with imf_option=1, for 2nd symmetric fission 9-May-2002 very minor bug in hauser-feshbach fixed. 10-April-2002 Got the program running on linux with LF95 (lehey-fujshu) this compiler is a bit more pickier than COMPAQ, so had to make a few changes to comform to its standard. 4-March-2002 forget to give time for fast-fission option - fixed 15-feb-2002 2) allowed the option to write out the residual excitation energy in a fragment at the point when GEMINI stopped decaying - this energy should go into gamma ray's. Use flag ex_flag to write these out. 2-feb-2002 1) modification to trig.f90 - better behavoured now 3-dec-2000 1) implimented full quantum treatment of fission angular distributions 3-Aug-2000 1)I have been making the code mode compatible with the F language and essential Lahey fortran- but similar subsets of fortran90/95. 2)Also for fast-fission, I allow the fragments to decay. 3)Added a new input parameter aden_type - give more flexibility for defining the temperature-dependent level-density parameter. 4)Compiled the code on an SGI machine and also with the LAHEY fortran90 compiler (had to remove some fortran95 features to do this) 26-Jan-2000 1) The R-matix expression I was using for 5He, 5Li and 8Be* line shapes had two much tail - This line shape is Eq 4 in Barker and Tracy Nucl. Phys. 38,33 (1962). I now had calculated these from the more correct expressions, eqs. 8,9 in the same reference. This make a big difference in the predicted yield of 8Be*. 2) I am getting impatient trying to get the multiplicities of of rare evaporated particles. I have therefore summed the probability of these emission at each step - and print out the multiplicities obtained from these probabilities at the end of the run. 19-Jan-2000 be8 array had a number wrong - be8 fist excited state gave wrong multiplicity 10-Jan-2000 fission delay option bug, - fixed 28 May-2001 found there was a error in cl_go - not sure when this occured, but its fixed now. introduced para%ssum again to speed things up when running fluctuations 24 Jan 2002 fixed a minor bug in write_mult ************************************************************** PROGRAM DESCRIPTION: ********************************************************************* CALCULATES THE DECAY OF A COMPOUND NUCLEUS VIA THE STATISTICAL MODEL ********************************************************************* This program follows the decay chain of a compound nucleus via sequential binary decays. The conditional barriers used in the statistical decay width for A < 100 were obtained from a two-spheroid finite-range program written by A.J. Sierk. For 100 < A < 190 the barriers are extrapolated from calculations of Sierk using more shape parameters. subroutine SADFIT performs a bi-cubic spline interpolation to obtain the conditional saddle point energy. ***** IMPORTANT ****** Use this program at your own peril, if you input angular momenta larger than where the fission barrier vanishes. For instance the subroutine of Sierks to give rotating ground state energies often gives rotation energies which decrease with angular momentum in this range. Sierk never intended his subroutines to be run at these angular momentum. Simular caution must be given concerning my extrapolation of asymmetric barriers. AUTHOR: Bob Charity Box 1134 Department of Chemistry Washington University St Louis Mo 63130 USA E-mail :- charity@wuchem.wustl.edu CREATION DATE: 1987 old version - 1999 fortran95 version ********************************************************************** INPUT: The program needs an input file of the form *.inp. example of the input file is as follows : ----------------------------------------------------------------------- Ni+zr liang ,title, 68 ,Z=Compound nucleus charge, 156 ,A=Compound nucleus mass, 113.94 ,Excn=excitation energy (MeV), 5.17 ,dexcn=excitation energy fluctuations (from target width) 0 ,lmin=minimum spin (hbar), 85 ,lmax=maximum spin (hbar), 10. ,R_num=*2l+1 events per l-wave, 0 ,N_weight=weighting factor, **************** output choices ************************ .FALSE. ,diagnostics, true=events typed out .false. .time_flag=.true. emission times written out on event file .false. .j_flag , fragment spins written on event file .FALSE. ,ex_flag, fragment excitation energies wriiten on event file .TRUE. ,I_angle=switch to calculate angles, .TRUE. ,true=quantum treatment of angles, false=semi-classical **************** fission and IMF parameters ************************ 1 ,imf_option,1=symmetric fission,2=all asymmetries,0=no imfs 0. ,t_delay=fission delay time in 10E-21 secndsseconds 0. ,sig_delay = first moment of delay as function of eta .TRUE. ,sharp_delay, gamma(t) = 0 (T) or gamma_0*t/t_delay(F) 3 ,Z_imf_min, minimum imf charge considered, 0.00 ,Kramers factor (only for IMF_option=1) .FALSE. ,Lestone - true =Lestone's summation over K included. 0. ,b_scale (default=1.00) **************** evaporation + gamma emission************************ .TRUE. ,tl_iwbc, true=tl's from IWBC model,False=sharp cut-off .FALSE. ,k_sum .FALSE. ,polarization 0 ,exotic_index 1.4 ,ratio 2 ,mass_option,0=liquid drop,1=experimental,2=shell fade out 50.0 ,E2 strength in Weiskofp units (default=50) 0.1 ,E1_strength 1.e-7 ,threshold *************** level density ***************************************** -6 ,aden_type, level density type 10. ,aden_0=level density constant (k=A/a) 1.02 ,a_scale ----------------------------------------------------------------------- I suggest you paste the above example in a file and modify it to your heart's desire. The Monte Carlo events are written in a event file called *.evt two ways to start GEMINI CONDITION - = new , the program starts a news calculation = continue, the program continues with a calculation which was interrupted by a computer crash or a premature program termination. In this mode, GEMINI reads the old event file for all input parameters and continues where it left off. No other input cards are required. ********************************************************************** The format of the input file is REACTION: - this can be the reaction title or whatever. The title is written as the first record in the output event file to help you identify it. Zcn: - Atomic number of compound nucleus Acn: - Atomic mass of compound nucleus Excn: - Excitation energy of compound nucleus. dexcn - fluctuation in the excitation (maybe useful for modeling a thick target experiment where there a range of excitation energies. Excitation are chosen randomly from excn-dexcn to excn+dexcn CN spin distribution The sharp cut-off approximation is assumed with. Lmin: - minimum compound nucleus spin Lmax: - maximum compound nucleus spin R_Num: - the number of Monte Carlo events per L-wave is (2*L+1)*R_Num N_weight: - = 0 no weighting. = 1 weighting scheme to produce lots of binary(fission) events. Gives good statistics for all mass splits. Weighting is turned off after a binary decay - so multi-fragment events may be rare. > 1 weighting kept on until N_weight binary decays (or fissions) occur. If you are interested in large fragment multiplicities, I suggest you use N_weight = 15 or maybe larger. If sig_delay >0, the code does not allow for weighting as of yet, and the code with stop. **************************************************************************** diagnostic = .TRUE. events are printed out as they are simulated useful for debugging purposes. .FALSE. nothing printed out on screen - option for major run. time_flag if .TRUE. the time at which a fragment was emitted is included in the information written to the event file. j_flag if .TRUE. the spin of the fragment is included in the information written to the event file. I_angle: if .TRUE. the angle and velocities of all the fragments are calculated if .FALSE. No angle and velocity information are calculated. **note when the angles are calculated, the spin vector of the compound nucleus is always assumed to point along the x-axis i.e. theta=90,phi=0. This allows one to quickly calculate the out-of-plane deflection. If a real experiment is to be simulated, then the phi angle can always be randomized in the sorting program. (or one can modify subroutine FUSION) NOTE angles will not be calculated at present with quantum=.true. quantum: if .FALSE. Semi-classical treatment of emission angles and velocities and spin orientations. (faster option) if .TRUE. quantum treatment of spin projections i.e. chooses the m values for each J. This option also give quantum angular distributions. I has not implemented this option for fission and IMF emission - *** note **** if you use this option the decay of the system will stop after the systems emits an IMF or fission. note both options are approximations and have their limitations 1)semi classical for evaporation - From the quantum mechanical numbers chosen in evaporation subroutines, classical vectors are determined. In the original version, the angle of the evaporated particle is always chosen perpendicular to the particle's orbital angular momentum vector - This however is only a good approximation for high angular momentum, and as most evaporated particles are emitted at low angular momentum, it produced too much emission in the reaction plane - For latest version, the angle of emission relative to the orbital angular momentum vector is obtained from the quantum distribution for m=l. 2)quantum treatment for evaporation - At each stage in the decay sequence, the "m" quantum numbers (angular momentum projection) are chosen for each angular momentum vector from the probabilities given by the square of Clebsch-Gordan coefficients. Once the "m" values are selected, the emission angle distribution is obtained from the square of the appropriate associated Legendre polynomial, the phi angles is always random. In the original version, the initial compound nucleus spin was assumed to be "m"=0 with quantization axis being the beam axis. This treatment is similar to that used by other evaporation codes such as EVAP. However, as all interference terms have been neglected in choosing the angle, there is no memory of the reaction plane. i.e. the phi angles are all random. As an alternative to this treatment, in the latest version the initial value of the compound nucleus spin projection is taken as "m"="j". Subsequently, all center of mass velocity vectors are rotated by theta=90 degrees and in phi by a single random angle. This produces almost identical angular distributions to the original treatment, but now the phi angles are no longer random but show some memory of the reaction plane, if the initial angular momentum is large. *************************************************************************** Fission and imf parameters imf_option: 0=no imf or fission decay 1=one fission mode only calculated - Bohr-Wheeler formalism this is identical to the codes PACE and CASCADE. In principle, this one fission mode shold include all asymmetries - however it does not lead to a prediction for the mass-asymmetry distribution - I therefore choose symmetric decay when choosing the fission fragments. If you are interested in fission mass distributions do noy use this option. It is appropriate for prediction cross sections. 2=all asymmetric divisions considered, Moretto's formalism *** note*** if you want sequential decay from the imf's and fission fragments, you must set quantum=.false. t_delay : if non zero, fission and Imf decay is inhibited for a time t_delay*EXP(-eta**2/(2*sig_delay**2))*1.e-21 seconds. if you select evaporation from a deformed - the system is considered deformed only during this period. (see ratio for more details) sig_delay :give eta=(A2-A1)/(A2+A1) dependence of delay time. (see line above) Z_imf_min :Choice of the dividing line between evaporation and and IMF emission. if Z=3, only Z=0,1,2 are treated in LIGHT_PARTICLE* routines with Hauser-Feshbach formalism and Z>=3 are treated in IMF with transition-state formalism. if Z=4, as well as Z=0,1,2 now Z=3 (i.e. 6Li+7Li with seven excited states) are treated in LIGHT_PARTICLE and Z>=4 are treated in IMF. This option is much slower. Kramer :scales down the fission decay width by a constant. if kramer=1, no scaling. Only works for imf_option=1. The Kramers factor should depend on the viscosity at the saddle point. Lestone : true - then the summation over K (projection of angular momentum on symmetry axis) is included when calculating the transition state decay width, (both fission and IMF emission) This was recently suggested by John Lestone. FOR SIMPLICITY, I have assumed the saddle-point moments of inertia are K independent, also I have assumed the ground state energy is also K independent. see Phys Rev C 59, 1540 (1999). : false - standard transitional state formalism USED b_scale : fission barrier = b_scale*sierk's fission barrier. only for imf_option=1. ******************************************************************************** Tl_iwbc ;switch to choose between sharp cut-off and incoming wave boundary condition (IWBC) model transmission coefficients .TRUE. = IWBC also the complete summation over orbital angular momentum and particle spin is performed .FALSE. = sharp cut-off, faster option Above summation not performed but decay widths are weighted by (2s+1) where s is particle spin. k_sum : if true, then the equilibrium distributions of K (projection of spin) and deformation are considered for calculating evaporation. See Phys Rev C61, 054614 (2000) for details The input value of ratio is ignored with this option. polarization : if true the L and M quantum numbers of all particles are written to the event file. Also transmission coefficients which depend on L.S (not averaged over) are used for protons and neutrons. exotic_index : selects the level of exotica of the evaporated fragments with Zratio>1 => oblate, ratio>1 =>prolate, ratio=0,1=> spherical The deformation is kept on only during the fission delay. Also in this peroid the yrast line used is that appropriate for the shape specified by ratio. after the fission delay, the transmission coefficients and yrast line appropriate for a "spherical system" are used. if you want to keep the deformation on for all evaporations, set imf_option=0 and t_delay>1000. I suggest you also use mass_option=0, as the "spherical" shell effects will be inappropraite for a deformed system. You will not get any fission or complex fragment emission. mass_option : =1 experimental masses are used to calculate all binding energies =0 Liquid drop masses of the parent and daughter nuclei and the experiment mass of the evaporated particle are used to calculate binding energies. (more appropriate at high excitation energies where shell effects are expected to wash out) =2 The shell corrections to the parent and daughter nuclei are assumed to fade out following the procedure of Reisdorf. E2_strength : the strength of E2 gamma rays in Weiskopf units default = 50 W.u. E+e1_strength : the strength of E1 gamma rays in Weiskopf units default = 10 W.u. threshold : meaningful only for tl_iwbc = .TRUE. Should be set to zero or a very small number. To save time when many exotic evaporation channels are considered, it is useful to turn off these channels as the excitation energy decreases to save CPU time. When the probability of a evaporation channel drops below the threshold value it is not calculated for the remainder of the evaporation cascade. I have used values of 1e-5 or less. ******************************************************************************* aden_type, aden_0 aden_0 is only used for aden_type=0,-9 Level-density constant parameters if aden_type = 0, then a constant level-density parameter given by little_a=A0/aden_0 is assumed if Aden_type=-1, Toke + Swiatecki level density parameters used, Nucl. Phys. A372 (1981) 141 if aden_type=-2, Ignatyuk et al level density (Yad Fiz 21 (1975) 1185 if aden_type=-3, Gottschalk+Ledergerber (Nucl. Phys. A278 (1997) 16 ****note ***** both (Toke + swiatecki), Ignatyuk, and (Gottschalk+Ledergerber) are mass asymmetry dependent and can affect the mass asymmetry distributions if aden_type=-4, Ormand et al temperature dependence of the level density parameter is used, However no surface area dependence, see Phys. Rev. C40 (1989) 1510. if aden_type =-6, lestones temperature dependent level density parameter, phys rev c52 (1995) 1118 Note on temperature dependent level density parameters. Temperature dependent level density parameters are are tricky things to deal with. One can define 3 temperature dependent little a's relating the temperature (T), thermal energy (u) and entropy (S) to each other. They are all different except when little a is constant. a(U,T) = U/T**2 a(S,T) = S/2/T a(S,U) = S**2/4/U In the statistical model, for calculating level densities, one wants to use a(S,U). However Shlomo and Natowitz calculations are for a(U,T). With GEMINI I generally use Lestone's (Phys. REV. C52 118 (1995) ) parameterization. He calculates a(S,U) , but it is a function of T which makes it difficult to implement as one has to iteratively solve for T first. After doing this a number of time, I found I could fit the final solutions and these are coded into GEMINI and used when aden_type = -6 is selected. aden_type = -9 - temperature-dependent little-a given by Fineman et al [Phys Rev C50, 1991 (1994)] a = A0/(8 + kappa*U/A0), where kappa = aden_0 a_scale : only used for imf_option=1. the level density parameter at the saddle-point configuration is scaled by this factor. if aden_type = 0, then a_scale = af/an, however other aden_type, the level density parameter may have a surface area dependence so af/an may already be different from unity and a_scale will further modify this value. Other parameter use may wish to change are Amax - The decay of products of A < or = to Amax is not considered. Exmax - the decay of products with excitation energy < or = to Exmax is not considered Z_shell - Shell and pairing corrections are included in the saddle point energies if the charge of the lighter product is < or = to Z_shell also for light-particle emission, the binding energies are calculated experimental mass for systems with A < or = Z_shell max_channels, max_frag - dimensions of arrays containing decay information for decay channel and all the emitted fragments. - at high excitation energies these may need to be increased. Program should tell you when you have trouble with these. These parameters need to be changed in the file GEMmod.f90 and the program needs to be recompiled and linked. The Monte Carlo events are written in the event file. To obtain charge distributions or whatever, one needs to sort this file. Program SORT is an example of a program to do this. Also Programs LOOK and FANCY allow you to look at event simulated event. These programs must be linked with SORTSUB. Program GEMINI Subroutines: Input : Reads the input file Event : Simulates the decay of one event DMP : Write a Monte Carlo event into the event file. Decay : Selects a decay channel via the Monte Carlo technique. decay_time_wise : Monte Carlo selection by time criteria decay_prob_wise Monte Carlo selection by probability criteria Decay_widths : supervises the calculation of decay widths Light_particle_* : calculates the light particle decay widths Imf : supervises the calculation of all complex fragment decay widths Complex_frag : calculates the decay width for a particular binary decay saddle : determines the saddle point energy gamma_ray : E1 and E2 decay widths fission_width : symmetric fission decay width non_statistical : decay of excited states of evaporated particles. at the moment only 6Li and 7Li decay treated. getmass : Supervises the calculation of a fragment mass Massfr : calculates the mass of a nucleus from the Yukawa-plus-exponential model Fusion : calculates the number of Monte Carlo events per L wave sadfit : bi-cubic spline interpolation of finite range model saddle-point energies jmax_sierk - gives maximum J for with symmetric fission barrier > 0. bfis_sierk - Sierks symmetric fission barrier yrast_sierk- Sierks yrast line M_inertia_sierk - moments of inertia of symmetric saddle point configuration. lpoly : used by *_Sierk subroutines yrast : gives the ground state energy of a nucleus as a function of spin yrastld : yrast in Liquid drop model sphere_imf : determine the spin and excitation energy of the fragments produced in a complex fragment decay and the emission angles and velocities of the fragments. sphere_lp: determine the emission angles and velocities of the light particle and residue. transform: performs a rotational transformation on spherical coordinates addvectors: adds two vector together. ran_j: selects the magnitude of the various collective modes via the MOnte Carlo method Gauss : calculates normal probability integral ********************************************************************* notes on weighting Complex fragment decays can be very improbable events. Thus in a normal Monte Carlo simulation one may have to run a large number events to achieve reasonable statistics for intermediate mass fragments. This will use a lot of CPU time and may cost a lot of money. To solve this problem, to some extent, weighting has been introduced in this program. If I_weight =1, the program includes weighting. At the start of an event, K_weight is set to 1. While K_weight =1 , the decay widths for complex fragment emission are increased relative to their calculated values. If in the decay chain, a complex fragment is emitted , then K_weight is set to 0 and the complex fragment decay widths are not weighted for the remainder of the decay chain. To take into account the increased yield of complex fragments , each event now has a weight. If a complex fragment is emitted the weight is decreased, whereas if no complex fragment is emitted, the weight is increased. The weight of an event is given as : Prob(statistical model)/Prob(simulation) Where Prob(statistical model) is the probability of that event, as calculated by the statistical model, and Prob(simulation) is the probability of GEMINI simulating that event. If we let GEMINI simulate events according to the statistical model probabilities (ie. I_weight=0), then the weight of each events is always 1, as is expected. If I_weight=1, Gemini selects the exist channels in the following manner: 1/ The complex fragment decay widths are multiplied by a constant (FACT) so that the first chance complex fragment emission probability is at least a value given by FACT0. (FACT0=0.1, but can be changed in gemini.inc) The same value of FACT is used for second, third, etc chance complex fragment emission. If the first chance complex fragment decay probability is already greater than "FACT0", then FACT is set to 1. Also after a complex fragment decay, FACT is set to 1 in simulating the decay of the primary fragments. Now the probability of gemini selecting a light particle decay channels is gamma(channels)/(gamma_total_light+Fact*gamma_IMF_total) where gamma(channel) is the decay width of the exit channel, gamma_total_light is the total decay width for light particle emission and gamma_IMF_total is the total decay width for complex fragment emission. Now the probability of choosing any complex fragment exist channel is of course FACT*gamma_IMF_total/(gamma_total_light+FACT*gamma_IMF_total. The actual channel is chosen so that all Charge-splits are equally likely. The mass split is then chosen according to the statistical model decay widths. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FREE ADVICE what does the message "entropy at sp > entropy at gs - 1" mean????????????? Let us first consider the original application of the transition state formulism - chemical reaction rates - There is a heat bath of constant temperature T and some molecule in this bath has an activation energy B (This is just the equilivant of the fission barrier) to disintigration. To derrive the transition state formulism, it is required that B/T >> 1. (1) You can look at at "reaction rate theory" Ha:nggi, Talkner, and Borkovec, Reviews of modern Physics V62 no2 (1990), to get an idea of why the condition is necessay (a bit detailed though). Anyway, this condition means the system is well bound as T is a typical value of the kinetic energy in the one dimensional "fission" mode and so the system does not have enough kinetic energy to surmount the barrier - it can only do so by exchanges of energy with the heat bath. Kramers treats the fission problem as a ramdon walk. Now for real nuclei, the temperature is not constant, it differs at the ground state and saddle-point. So I perfer a condition that does not involve the temperature directly. The condition should really be S(gs) - S(saddle point) >> 1. (2) where S is the entropy. Now S = (total energy - F)/T, where F is the free energy and the total energy is the same at the ground state and saddle point, hence equation [2] is equilivant to [F(saddle point) - F(ground state)]/ T(average) >> 1 where T(average) is an average temperture. This difference, F(saddle point) - F(ground state), is the definition of the temperature dependent fission barrier. note F(ground state) = E0 - an*T**2. and F(saddle point) = E0 + B - af*T**2. B is now the barrier at T=0. if "an" and "af" are not equal, then the fission barrier is temperature dependent. Generally it is assumed that af>an so the fission barrier decreases with temperature. Now in GEMINI, the code tells you when S(gs) - S(sp) = 1, but the real condition is >> not >. How much higher than 1 should you be, I am not sure, but it probability depends on how accurate you want to be. I note in the above reference, is it stated for chemical reaction rates, people generally only use the transition state formalism for B/T > 3. Another way of looking at this which might be more intutative is that the saddle point - is supposed to be a bottle neck through which the system has to pass to fission. The rate through the bottle neck determines the fission rate. - In a statistical theory, the population at any point is related to the level density which is proportional to EXP(S). If S(sp) > S(gs), then bottle neck is bigger than the bottle i.e. it can't be the rate limiting factor and so you can't apply the theory. Note all statistical model codes such as GEMINI, CASCADE, PACE, EVAP, MBII, ALICE, ALERT, etc, etc, use the transition state formalism - and are thus subject to condition [2]. It would not surprise me that other people are running these codes for systems where equation [2] is not valid. Most people do not appreciate all the approximation used in these formalisms and codes and run them as black boxes. ******************************************************************* more important notes _ LIMITATIONS OF GEMINI Extracted from a reply to a unhappy customer, who wondered why the latest version of GEMINI did not reproduce the predicted cross sections for Z=3,4,5 in the Nuclear Physics article. Let me first draw your attention to equation 10 in NP A483 (1988) p391. In the first version of GEMINI, I had a bug which gave incorrect values for M_exp in this equation. The modified saddle point energies were too large almost turning off Z=3,4,5 emission. In the present version of GEMINI, this bug has been eliminated and the yield are subsequently larger for Z=3,4,5. In fact Z=3 emission competes with Z=1,2 emission, which is why the evaporation residues have on average been shifted down in Z and are more spread out. Otherwise, the present version of GEMINI is consistent with the description in the above article. Now for your question as to what you can believe. First remember I am using Liquid drop (finite range corrected) barriers. Now also remember that the liquid drop model is never applied to very light nuclei when calculating nuclear binding energies, or at least the liquid drop model parameters are never fit to such nuclei. Thus for very asymmetric mass splits, i.e. when one of the nascent fragments is a very light nucleus, the predicted barriers are very questionable considering the model can't calculate the mass of this very light fragment anyway. Equation 10 in NP A483 was an attempt to correct this, the correction was only applied to mass splits where the lighter fragment was a Carbon or lighter. In my experience the corrected barriers are not good, but they are better than using the liquid drop values which would predict that Z=3 was the dominant decay mode in many cases. In fact the present code often predicts Z=3 emission is dominant at high angular momentum even with the correction. Thus to finally answer your question. Don't believe GEMINI for Z=3,4,5. For the larger Z's, the code does a reasonable job. It would not be unreasonable to turn off Z=3 as in my experience the code overestimates this. The residue cross sections would then be more reasonable. I am not sure turning off Z=4,5 would change any other cross sections. ***** in the newer version you have the option of using the Hauser Feshbach option to calculate the cross section for Z=3 and/or z=4 particles. This option is slower, but should give better results. ****************** old log of F77 version ******************************** C H A N G E L O G Date | Name | Description ----------------+-------+----------------------------------------------------- 14-5-89 | RJC | calculates the emission angles of all the fragments ----------------+-------+----------------------------------------------------- 6-10-89 | RJC | minor corrections to subroutine BARFIT and N_frag changed from LOGICAL*1 to INTEGER in GEMINI.INC ----------------+-------+----------------------------------------------------- 24-10-89 | RJC | modified the last check in subroutine DECAY_WIDTHS to be more in line with the PLI version. If I_channels>0, I also check to see if the total decay widths is greater than zero. This condition should never happen in reality - but due to underflows it can happen in the program - see notes embedded in the program. ----------------+-------+----------------------------------------------------- 5-12-89 | RJC | Major beautification - removed I_coll option - removed all double precision variables in BARFIT ----------------+-------+----------------------------------------------------- 8-6-90 | RJC | An extra weighting scheme has been introduced to produce more many-fragment events. The weighting flag I_weight is replaced by the variable N_weight. IF n_weight = 0 - no weighting = 1 old weighting (turned off after 1st binary decay) >1 weighting kept on until the N_weight binary decay. ----------------+-------+----------------------------------------------------- 20-8-90 | RJC | introduced a crude fission delay time into the program. All complex fragment channels are delayed by the same amount. The new input parameter is the delay time, if set to zero, the results are the same as the previous gemini version. Also added the IMPLICIT NONE statement to some routines and introduced a new format for the input file. ----------------+-------+----------------------------------------------------- 24-8-90 | RJC | Removed the few remaining labels from the code In doing so the dmp routine writes only once to the event file per event instead of twice. ***also the subroutine LIGHT_fragment is substantially modified. It no longer uses Sierk barriers to try an estimate a decrease in the Coulomb barrier with angular momentum. The inverse cross secnsections for n,p,a and their barriers come from the systematics of McMahan and Alexander Phys. Rev C21 (1980)1261 and Dostrovsky et al. Phy Rev 116 (1959)683. These are different to the original parameters indicated in the first GEMINI article. Warning, the new event file can longer be read by the old versions of SORT and LOOK. ----------------+-------+----------------------------------------------------- 30-8-90 | RJC | corrected two small errors of probably no consequence in subroutine barfit ----------------+-------+----------------------------------------------------- 31-1-91 | RJC | Major changes. The old version of subroutine LIGHT_PARTICLE which approximated the light particle decay widths by a classical integral over angular momentum - was not very good for low angular momentum (not surprising of course). There are now alternative subroutines which calculate the decay widths with varying degrees of accuracy. Light_particle = sharp-cut-off transmission coeff, but summation over angular momentum is performed correctly. This is what is described in original GEMINI paper. Light_particle_quick = sharp-cut-off transmission coeff where a summation over angular momentum is approximated - within 5% for all angular momentum for the cases I looked at. instead of summing over the variables l (particle orbital angular momentum) and s2 (residue spin). This routine sums over s2-l and then integrates over s2. Light_particle_very_slow = more realistic diffuse transmission coeff and correct summation over angular momentum. Allows for sub-barrier particle emission. The barriers and diffusenesses of the transmission coeff were approximated from the parabolic barrier approx, with the potential obtained from the real part of optical model parameterizations and then the curvatures and barriers were parameterized as simple functions of Z,A,l. ----------------+-------+----------------------------------------------------- 16-7-91 | RJC | 1) gamma ray decay branch added to terminate the particle when sub-barrier particle emission is allowed. 2) format of input file changed, the description of the parameters is listed after the input values. This saves one forgetting which line corresponds to which variable 3) The format of the event file has been changed, - both the Z and A of each fragment have been combined into one INTEGER*2 variable to save disk space 4) A subroutine fission_width is available to substitute for IMF if only symmetric fission is desired rather than all asymmetries. 5) new simplified treatment of light particle angular distribution in sphere_lp to go with new Light_particle subroutines ----------------+-------+----------------------------------------------------- 26-9-91 | RJC | the variable describing the properties of each fragment have all been included in a structured variable FRAG. the following are included in frag frag.Z atomic number frag.A atomic mass frag.Id (old I_decay) frag.Ex excitation energy frag.time decay time frag.spin = frag.spin.J = frag.j spin vector of fragment frag.spin.theta frag.spin.phi frag.vel = frag.vel.v velocity vector of fragment frag.vel.theta frag.vel.phi ----------------+-------+----------------------------------------------------- 21-11-91 | RJC | in choosing decay times assuming exponential decay we use LOG(RAN(seed)))/tau. Occasionally RAN(seed) is zero or small enough that this causes the program to crash. THis was prevented by adding 1.e-37 to RAN(seed) ----------------+-------+----------------------------------------------------- 12-1-93 | RJC | 1) Inclusion of a K xray decay mode. This is an atomic decay mode that does not complete against the nuclear decay modes, however as long the Z of the decaying system equals that of the compound system, an atomic K xray can be emitted. I assume that in the formation of the compound nucleus one k vacancy is created. This is definitely an over estimation - on average, in fusion reactions at low energies, values of less than a percent are expected, so the final k xray multiplicity must be multiplied by the K-vacancy production probability. No K xray emission is considered after the emission of a charged particles - such an emission can create K vacancies anyway which is too complicated to include and the anyway subsequent K xrays will have different energies. 2) I have made the selection of either "only symmetric fission" or "symmetric and all asymmetric mass splits" as an option to be read from the *.inp file. ----------------+-------+----------------------------------------------------- 20-1-93 | RJC | New VAX Fortran compiler complains of un-natural alignment in common blocks. Put all the INTEGER*2 variable at the end of the common block to get rid of these complaints. ----------------+-------+----------------------------------------------------- 7-4-93 | RJC | Gemini wasn't conserving energy. Some terms were forgotten in obtaining fragment kinetic energies in light_particle subroutines. This affected the velocity of the fragments - otherwise everything else Ok. Bug may be from 31-1-91 modification. ----------------+-------+----------------------------------------------------- 7-7-93 | RJC | It compiled on the our new DEC alpha without any modifications ------------------------------------------------------------------------------- 14-1-94 | RJC | 1) Made the level density constant k (or aden_0 in the code), where a=A/k, an input parameter, if k=-1 the Toke + Swiatecki NP A372 (1981) 141 is used, note this level density parameter is dependent on the mass split for fission-like decay. k=-2 Ignatyuk level density parameter is used k=-3 Gottschalk+ledergerber 2) included an input parameter to switch between IWBC and sharp cut-off transmission coeff. 3) The fission delay time is now dependent on mass asymmetry, delay=t_delay*EXP(-eta**2./(2.*sig_delay**2.) eta = A2-A1/(A2+A1) is mass asymmetry. Had to make extensive modifications to subroutine decay to achieve this. ----------------+-------+----------------------------------------------------- 6-4-94 | RJC | 1/ modified light_particle_very_slow to explicitly consider the spin of the emitted particle when integrating over l and residual spin. This option obtained with ssum=.true. Previously, the particle spin was ignored except for a degeneracy factor 2*spin+1. 2/ All light_particle evaporation subroutines can now evaporated 6Li and 7li in their ground and some excited states. These are obtained with mode set from 7 for 15 3/ Using newer subroutines of Sierk for ground state energies etc, also if angular momentum is greater than the value where the fission barrier vanishes, the yrast line is extrapolated up using a rigid momentum inertia obtained from Sierks MOMFIT routine. The Sierk subroutine gave strange results as it was never meant to be run for these angular momentum. 4/ a diagnostic flag is included so GEMINI prints out events as they are calculated. 5/Instead of stopping when a nucleus with fissility greater than that for 194Hg is encountered, the asymmetry dependence for 194Hg is assumed and the barriers are scaled so that the symmetric barrier is correct for the nucleus. THis is only useful if the nucleus is not too much heavier than 194Hg. A message is still printed out to warn you. 6/ As pointed out by Juan Cabrera, the A_cut parameter was never used in the code; IMF emission was allowed whenever possible. I modified the code so that now A_cut is used as in the original version and as described in the notes. ----------------+-------+----------------------------------------------------- 10-5-94 | RJC | 1/Using a new parameterization of the IWBC transmission coefficients. i.e. Tl = 1./(1.+exp(ee)) where ee = a1*ek * a2 + a3/ek + l**2.*(b1*ek + b2 + b3/ek) where ek - kinetic energy and l is orbital angular momentum, the a1,a2,a3,b1,b2,b3 constants were obtained from fitting. This parameterization has a better behaviour as ek->0 2/The Coulomb barriers for the sharp cut-off transmission coeff options are now obtained from the IWBC model where Tl = 0.5 to be more consistent with the above option. 3/Included option for quantum treatment of the spin projection for evaporation as opposed to the semiclassical treatment in sphere_lp. Obtained Clebsch Gordan coeff generator from Alex Brown at MSU. The subroutine quantum_lp could be extended to calculated the emission angles of the evaporated particles in the future. However, there is no quantum treatment coded for IMF emission. therefore only run quantum=.true. when IMF_option=0 ----------------+-------+----------------------------------------------------- 1-6-94 | RJC | 1)Include an option for fission delay, where the decay width increases linearly during the delay period reaching the equilibrium value. New input parameter - sharp_delay = FALSE for this option 2)Compiled with open VMS 6.1 for the first time, and got lots of warning messages about declared variables not being used. I removed these variables. ----------------+-------+----------------------------------------------------- 5-1-95 | RJC | Dave Bowman found a bug which can cause an array out of bound error in yrastld. Fixed. This subroutine is only rarely used so no big deal. ----------------+-------+----------------------------------------------------- 7-7-95 | RJC | summary of changes. 1. proper quantum mechanical treatment of angular distributions of evaporated particles is now included. use option QUANTUM = .TRUE. this is very slow If you run at high angular momentum you may need to increase some of the array sizes. i.e. increase l_max in quantum_ang_generator. If you never intend to use this option, you may want to decrease l_max to cut down on memory. 2. During fission delay, you can have evaporation from a deformed system (prolate). approximate transmission coeff for such a deformed system are available - I have commented them out in light_particle_very_slow as they really slow down the calculation. 3. There was an minor error in calculating the level density parameter from the Toke+swiatecki and other recipes, the wrong A was used. Also for aden=-4 , the temperature dependent level density parameter of Ormand is used ----------------+-------+----------------------------------------------------- 27-11-95 | RJC | 1/ a slower more correct subroutine for determination of the fission decay width was added. 2/ New input parameter - a Kramers factor which just scales down the fission decay width by a constant. 3/ Correction of a rare subscript out of range problem created during last modification. ----------------+-------+----------------------------------------------------- 9-12-95 | RJC | 1/ fixed bug created during last modification 2/ modified the K Xray code to allow K Xray emission after the evaporation of a light particle. If the system has a binary decay, K Xray emission is not followed after this time. In the event stream a K x-ray is identified as Z=0, A=-(Z of atoms). In the event file the minus sign gets lost and so A=256-(Z of atom) ----------------+-------+----------------------------------------------------- 13-12-95 | RJC | totally bored during an owl shift - spell checked the comments ----------------+-------+----------------------------------------------------- 16-2-96 | RJC | 1) For light_particle_very_slow, I changed the parameterization of the transmission coefficients, i.e. replaced the L term by an (L+1) term. I obtained slightly better fits to the IWBC transmission doing this. 2) Before quantum mechanical treatment of angular distributions was only allowed if all fission decay modes were turned off. i.e. IMF_option=0 Now the angular distributions for evaporated particles is calculated until a fission-like decay occurs in the event. After this no angles are calculated. This is useful if one wants only the angular distributions of light particles in coincidence with evaporation residues. ALSO increased the dimension of FCT in cl_go to allow higher l waves to be calculated 3) new options concerning the fragment masses. mass_option=0 : liquid drop masses are assumed as in all previous versions of GEMINI mass_option=1: experimental masses are assumed mass_option=2: experimental masses with a fading out of the shell correction with excitation energy is used. 6-11-96 | RJC | 1/ G_float option for quantum_ang_generator 2/ extended cha1 and cha2 and changed length of sad_array 16-1-97 | RJC | 1/ option to consider the L.S dependence of the transmission coefficients for protons and neutrons. This option is turned on if polarization=.true. This was included to study the polarization of protons emitted from a rotating compound nucleus of aligned spin. The L and M quantum number are written out for each particle when this option is selected 2/ Full treatment of deformation considered when the variable "ratio" is greater than unity, evaporation form a deformed compound system is simulated. The yrast lines is modified to account for the deformation and rotation energy of a prolate system with "ratio" being the ratio of major to minor axes. Also, depending of the level density option, the level density parameter may be modified to account for the increased surface area. Finally, the transmission coefficients are calculated following the procedure of Huizenga - averaging over the surface of the system. The system is assumed to remain deformed for a period "t_delay" and no fission competition is considered until after a spherical shape is assumed. Note the approximate deformed transmission coefficients mentioned earlier has been removed. 3/ The Hauser-Feshbach formalism can be used for the evaporation of all H, He, Li and Be isotopes including excited states with excitation energy less than 5 mev. Unstable states are allowed to sequentially decay via their known decay modes. There is a new variable "exotic_index" and the old variable Z_imf_min allow one to control which states you wish to include. Note, the Hauser-feshback formalism is much more time consuming than the transition state formalism. So unless you are interested in things like isotope ratios, some limitation of these states is suggested. Each state has been assigned an "exotic" number (see LIGHT_PARTICLE.INC), the larger the value, the more exotic is the decay. The variable exotic_index selects the level of exotica one wishes to include. Z_imf_min, selects the lowest Z value for which the transmission state formalism is applied. The Hauser_Feshback formalism is applied for all emissions with Z_imf_min-1, even if imf_option=1 and only symmetric fission or imf_option=0 and no decay are treated via the transition state formalism. 28-2-97 | RJC | fixed some bugs which very occasional cause the code to go into an infinite loop or crash. 7-4-97 | RJC | fixed a few more bugs which occasionally crashed the program 30-4-97 | RJC | 1) for mass_option=1,2 pairing effects are now assumed to be washed out 2) for quantum-treatment of angular distributions - the compound nucleus is now always assumed to the in an M=J state. Later all angles are rotated by 90 degrees so that the initial direction is at 90 degrees to the beam axis and reaction plane is randomly rotated. Previously we started with M=0, however this calculation produces no phi-correlation between emitted particles - i.e. no memory of the reaction plane. 3) introduced the new input variable a_scale to scale the level density parameter for the saddle point configuration. Beware, although a_scale can be used to set af/an, if aden_0 < 0, af/an is already not unity and a_scale further modifies it. However if aden_0 > 0 , then a_scale is af/an. 15-5-97 | RJC | 1) option tl_iwbc = .FALSE. was causing the program to crash - fixed this up, plus some other minor bugs. ************ this version now on FTP site for general distribution ******* 8-12-97 | RJC | fixed a bug with the collective enhancement option - gives realistic mass distributions. 17-12-97 | RJC | found for light systems, no p,d,t emission in huaser_Feshbach subroutine. the Fix, at "start just below the barrier" instead of ek_min = ek_min - 5 now start at ek_min = 0.8*ek_min 23-2-98 | RJC | improved the fit to the neutron transmission coefficients - not a major change Note for neutron with polarization turned on the optical model potential comed from Klepatski Phys Rev C34 2077 (1986) which unlike the stand Wilmore and Hodgson has a spin orbit term Note the fit for the transmission coefficient is bad for z<10, so don't use the polarization option here. 31-3-98 | RJC | 1) created new array para to make it easier to add new input parameters 2) added two new input parameters - dexcn and e2_strength and removed ssum and a_cut. 3) improved determination of deformed transmission coeff with quadratic interpolation instead of linear. 12-7-98 | RJC | 1)made subroutine transform more efficient, 2)made subroutine sphere_lp, semi classical angular distributions, more consistent with the quantum distributions 3) disabled the collective enhancement option 4) made the code more UNIX friendly 5) included option to calculate Lestone' transition state formalism which includes a summation over K 21-6-99 | RJC | 1) found and fixed a bug associated with the option shell fade out 2) made GEMINI Y2K compatible - i.e. the compiler no longer complains about Y2K potential problems which did not exist. 3) added some more Be and LI states one can evaporate. 7-6-99 | RJC | 1) the deformed option had been disabled- got it going again. 2) in previous version of gemini, the thermal excitation at saddle is divided between the two fragments - any energy dissipated in the saddle-to-scission motion is ignored. now for imf_option=1 (only), I add this extra energy to the excitation energy of the fragments. ************ this version now on FTP site for general distribution *******