The GetMech Input File
Explanation of the input file for getmech:
the following input file is included in the getmech.tar-ball. It has several entries. Here it is explained line for line.
Note, this input file is for a reaction that includes 4 intermediates plus the dark state. The reaction is followed by
27 time-points on a time scale from 1 ns to 4 s. The time-points are arranged liniarly on a logarithmic time-scale.
The crystal system (spacegroup) is hexagonal P63.
1 TIMEPOINTS, INTERMEDIATES, MASK MODE 1 RECOMMENDET
2 27 4 21
3 'extended_mask.pdb' 0.5 2 2.5 0
4 0
5 CELL CONSTANTS AND SYMMETRY OPERATORS:
6 66.9 66.9 40.8 90.0 90.0 120.0
7 6
8 X,Y,Z
9 -Y,X-Y,Z
10 Y-X,-X,Z
11 -X,-Y,1/2+Z
12 Y,Y-X,1/2+Z
13 X-Y,X,1/2+Z
14 TIME DEPENDENT MAPS, XFFT-FACTOR, LASER CORRECTION :
15 1E-9 '../MAPS_14/pyp14_1ns-dark_M30C_FS.map' 6055.5 1.00
2E-9 '../MAPS_14/pyp14_2ns-dark_M30C_FS.map' 5257.9 1.00
4E-9 '../MAPS_14/pyp14_4ns-dark_M30C_FS.map' 5003.3 1.00
8E-9 '../MAPS_14/pyp14_8ns-dark_M30C_FS.map' 4021.1 1.00
16E-9 '../MAPS_14/pyp14_16ns-dark_M30C_FS.map' 3591.4 1.00
32E-9 '../MAPS_15/pyp15a_32ns-dark_M30C_FS.map' 3594.0 1.00
64E-9 '../MAPS_14/pyp14_64ns-dark_M30C_FS.map' 3517.3 1.00
128E-9 '../MAPS_15/pyp15a_128ns-dark_M30C_FS.map' 3517.3 1.00
256E-9 '../MAPS_14/pyp14_256ns-dark_M30C_FS.map' 3484.9 1.00
512E-9 '../MAPS_15/pyp15a_512ns-dark_M30C_FS.map' 3484.1 1.00
1E-6 '../MAPS_14/pyp14_1us-dark_M30C_FS.map' 3238.5 1.00
2E-6 '../MAPS_15/pyp15a_2us-dark_M30C_FS.map' 3238.8 1.00
4E-6 '../MAPS_14/pyp14_4us-dark_M30C_FS.map' 3650.6 1.00
16E-6 '../MAPS_14/pyp14_16us-dark_M30C_FS.map' 3284.8 1.00
32E-6 '../MAPS_15/pyp15a_32us-dark_M30C_FS.map' 3253.0 1.00
64E-6 '../MAPS_14/pyp14_64us-dark_M30C_FS.map' 3219.8 1.00
256E-6 '../MAPS_14/pyp14_256us-dark_M30C_FS.map' 3351.5 1.00
512E-6 '../MAPS_15/pyp15a_512us-dark_M30C_FS.map' 3361.1 1.00
1E-3 '../MAPS_14/pyp14_1ms-dark_M30C_FS.map' 3410.0 1.00
4E-3 '../MAPS_14/pyp14_4ms-dark_M30C_FS.map' 3506.2 1.00
8E-3 '../MAPS_15/pyp15a_8ms-dark_M30C_FS.map' 3640.0 1.00
16E-3 '../MAPS_14/pyp14_16ms-dark_M30C_FS.map' 3714.2 1.00
32E-3 '../MAPS_15/pyp15a_32ms-dark_M30C_FS.map' 3831.1 1.00
128E-3 '../MAPS_15/pyp15a_128ms-dark_M30C_FS.map' 4015.2 1.00
512E-3 '../MAPS_15/pyp15a_512ms-dark_M30C_FS.map' 4907.9 1.00
2.0 '../MAPS_15/pyp15a_2s-dark_M30C_FS.map' 5370.1 1.00
41 4.0 '../MAPS_15/pyp15a_4s-dark_M30C_FS.map' 5266.0 1.00
42 TIME INDEPENDENT MAPS
43 'IMAPS/PIT_2PHY_FS.map' 861.620
44 'IMAPS/I0_2PHY_FS.map' 209.444
45 'IMAPS/pR_2PHY_FS.map' 400.469
46 'IMAPS/pB_2PHY_FS.map' 272.996
47 6
48 0.0362
49 2.5E8
50 2.8E7
51 30.4
52 0.639
53 11.4
54 'RECON/T15_recon'
55 'RECON/T15_resid'
line 1: ignored
line 2: 27 time points, 4 intermdiates, mask mode 21 mask mode 21 means: use pdb-file as mask and evolve the mask using a certain sigma value, you can use also small spheres to explicitly exclude some volume for analysis.
line 3: pdb-file around with the mask is calculated, margin in addition to Van-der-Waals Radii, 2 means that mask is evolved, 2.5 is the sigma value above and below grid points are allowed, 0 sigma is subtracted from the difference values.
line 4: 0, here no small spheres are used to exclude (all the evolved mask will be used).
line 5: ignored
line 6: cell constants
line 7: number of symmetry ops, here 6
line 8-13: symmetry ops in standard notation
line 14: ignored
line 15: time point, map file in quotes ‘ ‘, xfft-factor, Laser correction
the xfft-factor is needed to bring the difference maps in FSFOUR format back to the absolute scale the Laser correction can be used if the laser power varies from crystal to crystal.
line 41: last time point.
line 42: ignored
line 43-46: time-independent difference maps of the four intermediates
line 47: number of FIT-parameter (Nrates+1)
line 48: scale factor (1st fit parameter) = extent of reaction initiation
line 49-53: rate coefficients of the mechanism programmed into subroutine EVAL
line 54: file-name of calculated (reconstructed, mixed) time-dependent difference electron density maps
line 55: file-name of delta-delta-rho (difference-difference) maps = residual difference difference electron density between observed and calculated difference maps.
Programming of a New Mechanism into GetMech
Actually, the coding is relatively simple and only sounds more difficult than it actually is. You need to setup a scheme, that connects the intermediates, write down the coupled differential equations that account for this scheme and copy and paste this into the program.
Imagine you want to investigate the mechanism depicted below. PI0 is the first intermediate after the laser pulse. It decays to PR, PB and finally back the dark state PG. However, there are also other pathways possible characterized by rate coefficients k4 and k5. How to implement this into ‘getmech’?
C PI0 -> PR -> PB -> PG
C | k1 | k2 k3
C | \________>
C | k4
C |
C \________________>
C k5
Firstly, make a backup of your current ‘getmech’ version. Then open ‘getmech_v5.f’ and locate the subroutine EVAL1. It contains a kernel of a general, however very simple, coupled differential equation solver. The secret of this integration is that the small time period dt is dependent on the time-scale on which you are actually calculating. So, if you are on the ns time-scale, dt is some percentage (0.3%) of this time-scale (so dt is in the range of pico-seconds). If you progressed to the ms-time scale dt will be micro-seconds. This strongly avoids any rounding errors that could emerge if you calculate with pico-second sized dt into the second time-scale. In addition it makes this subroutine lightning fast. We will discuss here the kinetic mechanism that is shown above (it is also include in the distributed version). If you plan to implement a new mechanism, code it into EVAL1 and change EVAL1 to EVAL, since subroutine EVAL is used in the program (not EVAL1). Rename the other, older, EVAL subroutine to a name of your choice.
in EVAL1:
locate the line with EPS=1E-8
and TM=TM+DT
In between these two lines you just located you step with dt through the coupled differential equations. The array rates(N+1) contains the fit-parameter. rates(1) contains the overall scaling factor of the fit. It is also a fit parameter that is determined from the concentrations (the best-fit) of all intermediates to the X-ray data. rates(2) to rates(N+1) contains the N rate coefficients of the mechanism. If you have a mechanisms with N different from the one depicted here, you only need to change the rate assignment and replace the differential equations below with new ones.
PI0 decreases in favor of PR with k1 and of PG with k5
dPI0 = ( -PI0 * k1 – PI0*k5 ) dt
PR gains from PI0 and loses to PB with k2 and loses to PG with k4
dPR = (+PI0 * k1 – PR*k2 – PR*k4) dt
PB gains from PR and loses to PG with k3
dPB = (+ PR*k2 – PB*k3 ) dt
pG gains from all intermediates:
dPG = (+PI0*k5 + PR*k4 + PB*k3) dt
PI0=PI0+dPI0
PR=PR+dPR
PB=PB + dPB
PG=PG+dPG
Care is taken in the program that you won’t get negative concentrations.
The calculated concentrations for all intermediates are then stored in the array TA. You need to locate this array at the end of EVAL1. If you have NTP time-points, the concentrations of the first intermediate need to be stored in TA(1 to NTP), the concentrations of the second intermediate will be stored in TA(1+NTP … 2*NTP) and so on. If you have more intermediates than in this example for example 4 intermediates, you need to store the concentrations of the 4th intermediate in TA(III+3*NTP). III is the actual time-point. If you have less intermediates you need to modifiy TA accordinly. Concentrations in array TA are used to mix the calculated time-independent differerence electron density maps. EVAL is also used to write out the concentration to standard output.
IF (IWRITE.GT.0) THEN
WRITE(*,500) TM,PI0,PR,PB,PG
500 FORMAT(G15.3,4F10.5)
make sure to change the names of the intermediates as well as the FORMAT statement when you implement a new mechanism.
This is all you need to take care of. Everything else should be automatic and/or specified by the input-file.
EVAL is then called by a routine called DDRHO that generates time-dependent difference electron densities from input time-independent difference electron densities of the intermediates and calculates the difference to the observed difference electron densities. Of course, the goal is to minimize this difference by changing the fit-parameters namely the scale and the (here 5) rate coefficients. This is done by the main fit-driver routine called FIT. If everything goes well, ‘getmech’ performs the fit.