diff --git a/trimsp/src/Changelog.txt b/trimsp/src/Changelog.txt new file mode 100644 index 0000000..3d13fbb --- /dev/null +++ b/trimsp/src/Changelog.txt @@ -0,0 +1,169 @@ +Sep 1998: to create an executable for PC running under WIN95 + using the DIGITAL VISUAL FORTRAN compiler (ed. Aug. 97) + Insertion of !DEC$REAL:8 (all REAL are REAL*8) + Conversion of function to double precession + Conversion of function to integer*4 + Conversion of random number generator from DRAND48() + to DBLE(RAN(ISEED)) + UNIT 17 (data to tape) disabled + UNIT 11 (input data) changed to file name EINGABE1.inp + UNIT 21 (output.data) changed to file name AUSGABE1.out + Insertion of UNIT 22 (only range data), file name AUSGABE1.rge +Dec 1998: Introduction of gaussian distributed projectile energies + new variable RI2 = random number initializer for gaussian energy distribution + new variable ISEED2 = random number for gaussian energy distribution + new variable Esig = sigma of the energy distribution + if Esig=0. then fixed projectile energy + new variable Epar = projectile energy + new subroutine ENERGGAUSS = for the gaussian distribution + new variables in subroutine p1,p2,p3, necessary for calculation of + the gaussian energy distribution +Mar 1999: LOGICAL FUNCTION EQUAL inserted. This function is used for comparision + of REAL*8 variables with 0. + OUTPUT to files AUSGABE1: + bug fixed for the output of the CH(i), bug was also present + in the original and the older version of TRIMSPP*C + if OUTPUT file exists then program asks for new OUTPUT filename + all variables are now defined as REAL*8 or INTEGER*4 + for variable conversion from REAL to INTEGER or INTEGER to REAL + conversion function inserted (IDINT, DFLOTJ) + use of correct MIN or MAX functions + data initialization for all variables (like in BHam) but without + DATA ier/102*0/ give an error in calculation of ions stopped in layer 2 - 3 + Unit 99 for file ausgabe1.err inserted, if and IF... (see below is true + then a message is included into this file + IF's inserted - can be found by using find WRITE(99, + DABS inserted for DLOG and DLOG10 arguments +Mai 1999: Version h + for TRIMSP simulations running in batch mode (i.e. to calculate + a serie of different range profiles using the same target an output + file FOR33 is created. In this file the following parameters are + written: + E0,Esig,NH,IIM,IB,IT,IRPL(1),IRPL(2),IRPL(3),fix0,sigmax + the FILE open statement is inserted after label 1502 + new variable column(100)*104 and colcount inserted. +Mai 1999: Version i + tried to fix bug for variable C2 (in Do 65 loop). + Sometimes this variable becames smaller then 10-10 (an underflow occurs, + the variable C2 is NaN, variables S2 and DENS, which are directly calculated + from C2 are also NaN. If variable C2 is less then 10-10 then C2 is 10-10. + S2 and DENS are calculated with the resetted value of C2. + Range file (output to unit 22) slightly changed. Now the midth of a channel + is calculated and only the number of particles is written to the file + +Jun 1999: Version j1 + changed calculation of Firsoc screening length according to the IRCU 49 report. + the Andersen-Ziegler tables (file stopping.dat) are revised,i.e. the coefficients + from the ICRU 49 report are included for all elements (file stopicru.dat) + + PROGRAM DATMAK2c: this program creates the input files for TRIMSPP4*. + with this program one can change the energy, + the energy distribution, the number of muons, and the + layer thickness easily. This is for making different + simulations (i.e. energy scan, sigma scan, particle + number scan, layer thickness scan). A corresponding + batch file is created with this program too. +Nov 1999: Version j1a + UNIT33: now Etrans, sigmaEtrans, Eback, sigmaEback included + Bug fixed in line 2084, location of 'write to unit33' changed. +Dec 1999: Version k + inclusion of variables + tryE :how often a random energy is calculated by the random generator + for large E0 this number should be the same as the number of projectiles nh + negE :how often a negative energy is calculated by the random generator + tryE - negE should be nh + both variables are type integer, both variables included in the output file + and in the fort33.file + inclusion of alphasig: one dimensional distribution of the angle of incidence + if alphasig .NE. 0 THEN a gaussian distribution of alpha will be calculated + in the subroutine ALPHAGAUSS + if the calculated angle is < 0 then the absolute value is taken because of + the options alpha = -1 or alpha = -2 (see file TRVMC95-v4k.txt) + inclusion of RI3: random seed for the calculation of the gaussian distribution + of alpha. + Header line for file for33 included. +Jan 2000: no new version but + included energy scaling of particle reflection coefficients after Thomas et al. + NIM B69 (1992) 427 + the calculated particle reflection coefficients prcoeff are written to the fort.33 file + prc is calculated if 1st layer consists only one element ! + included constants prc(1) - prc(6) + Thomas-Fermi reduced energy: epsilon + output of E0 and Esig now in keV +Apr 2000: BUG fixed, Var CHM2 removed, i.e. stopping power for energies below 10 keV + was a little bit overestimated. now version TrimSpP4L + Prg Datmak4L creates input files. +Jun 2000: Source code enlarged to simulate implantation profiles in up to 7 layers. + Each layer can have up to 5 different elements. + BUG fixed in line 1470,i.e. in calculation of scattering angle for the + Moliere potential the variable rrr1 was not handled correct. The length of + line 1470 was too long, RRR1 was read as RRR which has the value 0. + Due to the fact, that all calculation based on the other interaction potentials + (Krypton Carbon and ZBL) give more or less the same inplantation profiles, this + bug has a minor influence to the implantation profile. + Variables DateTime(8) = Integer array and + Real_Clock(3) = Character array included for date and time output to file + Changes for output file to UNIT 33: + Size of character variable column enlarged from 214 to 246 + Program datmak7a creates input files for this TrimSp release. + true calculation, i.e. indendent of bin width, of particles stopped + in layers - see UpTiefe,LowTiefe,number_in_layer(7) + INTEGER check_layer_flag : IF 1 then calculated implantation profile + in the 100 depth intervals agrees with + number of particles stopped in the different layers + IF 0 not, message will be written in the range file + UNIT 21 and UNIT22 +Nov 2000: Tanya Riseman + Starts porting program to Open VMS and DEC UNIX. + Compile options: f90/list/warn/ext + Minor problems with continuation characters in wrong column and + with a few variables undeclared, mostly functions. + Make code fit on 72 columns, because -extend_source + does not appear to work on DEC Unix. + Start changing strings in format and write statments so that + strings don't straddle continuation character in + column 6. Straddling can be non-portable! + Comment out !DEC$REAL:8 because all + reals are already declared REAL*8. Add implicit none. + +Jun 2002: Thomas Prokscha PSI + Stopped porting to f90, use f77 only. + replaced all non-standard f77 function by standard functions. + use ranlux random number generator from the CERN library (libmathlib.a must + be installed) to get rid of machine specific random number generators. + Add pre-compiler instructions for making different output for the .rge files + on Windows and Unix. + +Oct 2002: Thomas Prokscha PSI + corrected error in the calculation of the Thomas-Fermi reduced energy: + it was 1/(Z1 Z2 * sqrt( Z1**2/3 + Z2**2/3)) + it must be: + 1/(Z1 Z2 * sqrt( Z1**(2/3) + Z2**(2/3))) + +Aug 2009: Zaher Salman PSI + Changed input file reading to non-formatted. This still works + with old input files, but makes it much easier to create new + input files without the hassle of formatting. + Included the source of ranlux for random numbers generation + into trimsp7l source. No need for cern libraries to be installed. + +Feb 2013: Zaher Salman PSI + Cleaned up the code. + When possible remove loop numbres and use do-enddo instead. + Using proper fortran line indentation (emacs style). + Created TimeStamp subroutine to generate time stamp and removed + original code from main (two places) + +Mar 2013: Zaher Salman PSI + Started implementation of a more flexible input file format, with + flexible number of layer, not limited to max 7. + +Apr 2013: Zaher Salman PSI + Implantation into structures up to 100 layers is possible + Depth profiles with up to 500 points are now possible + + Note: The .out file contains only the first 7 layers. This is due + Fortran formatting limitation and is unlikely to changes unless the + layout of the output files is changed. However, for most purposes + this does not affect the functionality of the simulation. + diff --git a/trimsp/src/trimspNL.F b/trimsp/src/trimspNL.F index f91db32..1e1bf8a 100644 --- a/trimsp/src/trimspNL.F +++ b/trimsp/src/trimspNL.F @@ -31,168 +31,6 @@ C ENERGY DISTRIBUTIONS IN SPECIFIC DIRECTIONS C ENERGY DISTRIBUTIONS IN 100 STEPS UP TO INCIDENT ENERGY C ANGULAR DISTRIBUTIONS IN STEPS OF COS = 0.05 (CONST. SOLID ANGLE) C -C MAJOR CHANGES: -C -C Sep 1998: to create an executable for PC running under WIN95 -C using the DIGITAL VISUAL FORTRAN compiler (ed. Aug. 97) -C Insertion of !DEC$REAL:8 (all REAL are REAL*8) -C Conversion of function to double precession -C Conversion of function to integer*4 -C Conversion of random number generator from DRAND48() -C to DBLE(RAN(ISEED)) -C UNIT 17 (data to tape) disabled -C UNIT 11 (input data) changed to file name EINGABE1.inp -C UNIT 21 (output.data) changed to file name AUSGABE1.out -C Insertion of UNIT 22 (only range data), file name AUSGABE1.rge -C Dec 1998: Introduction of gaussian distributed projectile energies -C new variable RI2 = random number initializer for gaussian energy distribution -C new variable ISEED2 = random number for gaussian energy distribution -C new variable Esig = sigma of the energy distribution -C if Esig=0. then fixed projectile energy -C new variable Epar = projectile energy -C new subroutine ENERGGAUSS = for the gaussian distribution -C new variables in subroutine p1,p2,p3, necessary for calculation of -C the gaussian energy distribution -C Mar 1999: LOGICAL FUNCTION EQUAL inserted. This function is used for comparision -C of REAL*8 variables with 0. -C OUTPUT to files AUSGABE1: -C bug fixed for the output of the CH(i), bug was also present -C in the original and the older version of TRIMSPP*C -C if OUTPUT file exists then program asks for new OUTPUT filename -C all variables are now defined as REAL*8 or INTEGER*4 -C for variable conversion from REAL to INTEGER or INTEGER to REAL -C conversion function inserted (IDINT, DFLOTJ) -C use of correct MIN or MAX functions -C data initialization for all variables (like in BHam) but without -C DATA ier/102*0/ give an error in calculation of ions stopped in layer 2 - 3 -C Unit 99 for file ausgabe1.err inserted, if and IF... (see below is true -C then a message is included into this file -C IF's inserted - can be found by using find WRITE(99, -C DABS inserted for DLOG and DLOG10 arguments -C Mai 1999: Version h -C for TRIMSP simulations running in batch mode (i.e. to calculate -C a serie of different range profiles using the same target an output -C file FOR33 is created. In this file the following parameters are -C written: -C E0,Esig,NH,IIM,IB,IT,IRPL(1),IRPL(2),IRPL(3),fix0,sigmax -C the FILE open statement is inserted after label 1502 -C new variable column(100)*104 and colcount inserted. -C Mai 1999: Version i -C tried to fix bug for variable C2 (in Do 65 loop). -C Sometimes this variable becames smaller then 10-10 (an underflow occurs, -C the variable C2 is NaN, variables S2 and DENS, which are directly calculated -C from C2 are also NaN. If variable C2 is less then 10-10 then C2 is 10-10. -C S2 and DENS are calculated with the resetted value of C2. -C Range file (output to unit 22) slightly changed. Now the midth of a channel -C is calculated and only the number of particles is written to the file -C -C -C Jun 1999: Version j1 -C changed calculation of Firsoc screening length according to the IRCU 49 report. -C the Andersen-Ziegler tables (file stopping.dat) are revised,i.e. the coefficients -C from the ICRU 49 report are included for all elements (file stopicru.dat) -C -C PROGRAM DATMAK2c: this program creates the input files for TRIMSPP4*. -C with this program one can change the energy, -C the energy distribution, the number of muons, and the -C layer thickness easily. This is for making different -C simulations (i.e. energy scan, sigma scan, particle -C number scan, layer thickness scan). A corresponding -C batch file is created with this program too. -C Nov 1999: Version j1a -C UNIT33: now Etrans, sigmaEtrans, Eback, sigmaEback included -C Bug fixed in line 2084, location of 'write to unit33' changed. -C Dec 1999: Version k -C inclusion of variables -C tryE :how often a random energy is calculated by the random generator -C for large E0 this number should be the same as the number of projectiles nh -C negE :how often a negative energy is calculated by the random generator -C tryE - negE should be nh -C both variables are type integer, both variables included in the output file -C and in the fort33.file -C inclusion of alphasig: one dimensional distribution of the angle of incidence -C if alphasig .NE. 0 THEN a gaussian distribution of alpha will be calculated -C in the subroutine ALPHAGAUSS -C if the calculated angle is < 0 then the absolute value is taken because of -C the options alpha = -1 or alpha = -2 (see file TRVMC95-v4k.txt) -C inclusion of RI3: random seed for the calculation of the gaussian distribution -C of alpha. -C Header line for file for33 included. -C Jan 2000: no new version but -C included energy scaling of particle reflection coefficients after Thomas et al. -C NIM B69 (1992) 427 -C the calculated particle reflection coefficients prcoeff are written to the fort.33 file -C prc is calculated if 1st layer consists only one element ! -C included constants prc(1) - prc(6) -C Thomas-Fermi reduced energy: epsilon -C output of E0 and Esig now in keV -C Apr 2000: BUG fixed, Var CHM2 removed, i.e. stopping power for energies below 10 keV -C was a little bit overestimated. now version TrimSpP4L -C Prg Datmak4L creates input files. -C Jun 2000: Source code enlarged to simulate implantation profiles in up to 7 layers. -C Each layer can have up to 5 different elements. -C BUG fixed in line 1470,i.e. in calculation of scattering angle for the -C Moliere potential the variable rrr1 was not handled correct. The length of -C line 1470 was too long, RRR1 was read as RRR which has the value 0. -C Due to the fact, that all calculation based on the other interaction potentials -C (Krypton Carbon and ZBL) give more or less the same inplantation profiles, this -C bug has a minor influence to the implantation profile. -C Variables DateTime(8) = Integer array and -C Real_Clock(3) = Character array included for date and time output to file -C Changes for output file to UNIT 33: -C Size of character variable column enlarged from 214 to 246 -C Program datmak7a creates input files for this TrimSp release. -C true calculation, i.e. indendent of bin width, of particles stopped -C in layers - see UpTiefe,LowTiefe,number_in_layer(7) -C INTEGER check_layer_flag : IF 1 then calculated implantation profile -C in the 100 depth intervals agrees with -C number of particles stopped in the different layers -C IF 0 not, message will be written in the range file -C UNIT 21 and UNIT22 -c Nov 2000: Tanya Riseman -c Starts porting program to Open VMS and DEC UNIX. -c Compile options: f90/list/warn/ext -c Minor problems with continuation characters in wrong column and -c with a few variables undeclared, mostly functions. -c Make code fit on 72 columns, because -extend_source -c does not appear to work on DEC Unix. -c Start changing strings in format and write statments so that -c strings don't straddle continuation character in -c column 6. Straddling can be non-portable! -c Comment out !DEC$REAL:8 because all -c reals are already declared REAL*8. Add implicit none. -c -c Jun 2002: Thomas Prokscha PSI -c Stopped porting to f90, use f77 only. -c replaced all non-standard f77 function by standard functions. -c use ranlux random number generator from the CERN library (libmathlib.a must -c be installed) to get rid of machine specific random number generators. -c Add pre-compiler instructions for making different output for the .rge files -c on Windows and Unix. -c -c Oct 2002: Thomas Prokscha PSI -c corrected error in the calculation of the Thomas-Fermi reduced energy: -c it was 1/(Z1 Z2 * sqrt( Z1**2/3 + Z2**2/3)) -c it must be: -c 1/(Z1 Z2 * sqrt( Z1**(2/3) + Z2**(2/3))) -c -c Aug 2009: Zaher Salman PSI -c Changed input file reading to non-formatted. This still works -c with old input files, but makes it much easier to create new -c input files without the hassle of formatting. -c Included the source of ranlux for random numbers generation -c into trimsp7l source. No need for cern libraries to be installed. -c -c Feb 2013: Zaher Salman PSI -c Cleaned up the code. -c When possible remove loop numbres and use do-enddo instead. -c Using proper fortran line indentation (emacs style). -c Created TimeStamp subroutine to generate time stamp and removed -c original code from main (two places) -c -c Mar 2013: Zaher Salman PSI -c Started implementation of a more flexible input file format, with -c flexible number of layer, not limited to max 7. c------------------------------------------- c check OS c @@ -204,11 +42,14 @@ c IMPLICIT NONE C These parameters are related to the maximum number of layers MAXNL -C and define the number of points in the depth distribution MAXD +C and maximum number of points in the depth distribution MAXD INTEGER MAXD,MAXD1,MAXD2,MAXD5,MAXDNL5,MAXD2p2,MAXD2meagb & ,MAXD2meab INTEGER MAXNL,MAXNL5,MAXNLp25,MAXNL5p2,MAXNLm15 - PARAMETER (MAXD=200) +C This is the only point where the number of layers and depth +C profile are changed. All other parameters shouold be changed +C accordingly. + PARAMETER (MAXD=500) PARAMETER (MAXNL=100) PARAMETER (MAXD1=MAXD+1) PARAMETER (MAXD2=MAXD+2) @@ -637,7 +478,7 @@ C open statement for output files, removed from line 2449 ff to here OPEN(UNIT=21,FILE=outnam,STATUS='new') 6001 OPEN(UNIT=22,FILE=rgenam,STATUS='new') WRITE(21,1000) - 1000 FORMAT(1H1/,6X,'* PROGRAM TRVMC95 - V TrimSP7L 17.Oct.02 TP *') + 1000 FORMAT(1H1/,6X,'* TrimSPNL 02.Apr.13 *') C Get simulation start time CALL TimeStamp(day_start,month_start,year_start,hour_start @@ -2825,7 +2666,14 @@ C #endif IF(YH.LT.1.D0) GO TO 603 DO I=0,MAXD1 +C This normalization of the implantation profile is wrong if the +C simulation does not account for all implanted projectiles, +C e.g. when the number of points in the profile is not enough to +C reach the maximum implantation depth. RIRP(I) = DBLE(IRP(I))/YH +C To resolve this issue we need to check whether we reached maximum +C implantation or not. If not rerun calculation with larger step +C size! ENDDO 603 D1=0. D2=CW