%%bash

# Align bound DNA
echo 6 0 | gmx trjconv -f inputs/F1_complex_DNA.xtc -s inputs/dna.tpr -n inputs/dna.ndx -o complex_dna_aligned.xtc -fit rot+trans

# Align free DNA
echo 6 0 | gmx trjconv -f inputs/F1_free_DNA.xtc -s inputs/dna.tpr -n inputs/dna.ndx -o free_dna_aligned.xtc -fit rot+trans

"""
Explanation: How to calculate Global Elasticity?
To calculate global elasticity, dnaMD globalElasticity command can be used. However, it takes HDF5 file as input.
Following steps can be performed to generate HDF5 files. The tutorial file can be downloaded here. We will prepare HDF5
file for both free and bound DNA.
Calculate stretching twisting and bending motions
Both free and bound DNA is superimposed on to the same DNA structure. Careful that bending calculation is fitting dependent.
Therefore, at first we aligned both free and bound DNA to a common DNA structure as follows
End of explanation
"""


%%bash

# For free DNA
echo 2 | do_x3dna -f free_dna_aligned.xtc -s inputs/dna.tpr -n inputs/dna.ndx -ref -noavg -nofit -name free
mv *_free.dat outputs/.

# For bound DNA
echo 2 | do_x3dna -f complex_dna_aligned.xtc -s inputs/dna.tpr -n inputs/dna.ndx -ref -noavg -nofit -name bound
mv *_bound.dat outputs/.

"""
Explanation: 1. Run do_x3dna on DNA trajectory, -nofit is used because DNA is already superimposed to a common DNA structure using trjconv.
End of explanation
"""


%%bash

# For free DNA
dnaMD saveAsH5 -tbp 27 -i outputs/L-BPS_free.dat,outputs/L-BPH_free.dat,outputs/HelAxis_free.dat -o free_dna.h5
dnaMD axisCurv -tbp 27 -bs 2 -be 25 -ctan -scp 100 -s 1000 -cta 30 -io free_dna.h5 -ap free_dna_axis.pdb

# For bound DNA
dnaMD saveAsH5 -tbp 27 -i outputs/L-BPS_bound.dat,outputs/L-BPH_bound.dat,outputs/HelAxis_bound.dat -o bound_dna.h5
dnaMD axisCurv -tbp 27 -bs 2 -be 25 -ctan -scp 100 -s 1000 -cta 30 -io bound_dna.h5 -ap bound_dna_axis.pdb

"""
Explanation: 2. Run dnaMD to extract the parameters from do_x3dna output files and save as HDF5 file. Also calculate global axis,
curvature and tangents.
End of explanation
"""


%%bash

dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -paxis X -o elastic_modulus_BST.csv

"""
Explanation: Now, we have HDF5 files of both free and bounds DNA. It can be used for the calculation of elastic properties. These files
can be used with either dnaMD Python module or dnaMD globalElasticity.
Bending Stretching Twist Modulus
Following command calculate Bending Stretching Twist modulus matrix. Output matrix will be stored in csv file.
Elastic modulus matrix is printed as output and average values of contour length and cumulative twist angle
is also printed.
End of explanation
"""


%%bash

dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype ST -o elastic_modulus_ST.csv

"""
Explanation: The above modulus matrix is in this form:
$$\text{modulus matrix} =
\begin{bmatrix}
M_{Bx}       & M_{Bx,By} & M_{Bx,S} & M_{Bx,T} \
M_{Bx,By}    & M_{By}    & M_{By,S} & M_{By,T} \
M_{Bx,S}     & M_{By,S}  & M_{S}    & M_{S,T} \
M_{Bx,T}     & M_{Bx,T}  & M_{S,T}  & M_{T}
\end{bmatrix}
$$
Where:

$M_{Bx}$ - Bending-1 stiffness in one plane
$M_{By}$ - Bending-2 stiffness in another orthogonal plane
$M_{S}$ - Stretch Modulus
$M_{T}$ - Twist rigidity
$M_{Bx,By}$ - Bending-1 and Bending-2 coupling
$M_{By,S}$ - Bending-2 and stretching coupling
$M_{S,T}$ - Stretching Twsiting coupling
$M_{Bx,S}$ - Bending-1 Stretching coupling
$M_{By,T}$ - Bending-2 Twisting coupling
$M_{Bx,T}$ - Bending-1 Twisting coupling

Stretching Twist Modulus
Following command calculate Bending Stretching Twist modulus matrix. Output matrix will be stored in csv file.
Elastic modulus matrix is printed as output and average values of contour length and cumulative twist angle
is also printed.
End of explanation
"""


%%bash

dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -paxis X -fgap 100 -em block -ot modulus_time_BST.csv

# Print first  and last 3 line of output file
echo "====================================="
echo "Elastic modulus as a function of time"
echo "====================================="
head -4 modulus_time_BST.csv
printf ".\n.\n.\n"
tail -3 modulus_time_BST.csv

"""
Explanation: The above modulus matrix is in this form:
$$\text{modulus matrix} =
\begin{bmatrix}
M_{S}    & M_{S,T} \
M_{S,T}  & M_{T}
\end{bmatrix}
$$
Where:

$M_{S}$ - Stretch Modulus
$M_{T}$ - Twist rigidity
$M_{S,T}$ - Stretching Twsiting coupling

Convergence in Modulus
Same command can be used to calculate elasticity as a function of time with option 
-ot/--output-time and save it in csv format file. This result can beused to check
their convergence. 

If this option is used, -fgap/--frame-gap is an essential option.
This options also gives final value and error as the output on display.
The output file is in csv format and can be opened as spreadsheet.

NOTE: Elastic properties cannot be calculated using a single frame because 
fluctuations are required. Therefore, here time means trajectory between zero 
time to given time.
End of explanation
"""


%%bash

dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype ST -paxis X -fgap 100 -gt "gmx analyze" -em "block" -ot modulus_time_ST.csv

# Print first  and last 3 line of output file
echo "====================================="
echo "Elastic modulus as a function of time"
echo "====================================="
head -4 modulus_time_ST.csv
printf ".\n.\n.\n"
tail -3 modulus_time_ST.csv

"""
Explanation: Some of the plots from above data can be found here

Same as above but only for stretching and twisting motions.
End of explanation
"""


%%bash

dnaMD globalEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -et "all" -paxis X -gt "gmx analyze" -em "block" -o energy_all_BST.csv

# Print first  and last 3 line of output file
echo "==========================================================="
echo "Deformation free energy of bound DNA as a function of time"
echo "==========================================================="
head -4 energy_all_BST.csv
printf ".\n.\n.\n"
tail -3 energy_all_BST.csv

"""
Explanation: Global deformation free energy
To caluclate global deformation enrgy, dnaMD globalEnergy can be used. At first, elastic matrix from reference
DNA (most often free or unbound DNA) is calculated and subsequently this matrix is used to calculate deformation free
energy of probe DNA (most often bound DNA).
The deformation free energy is calculated using elastic matrix as follows
$$G = \frac{1}{2L_0}\mathbf{xKx^T}$$
$$\mathbf{x} =  \begin{bmatrix}
                (\theta^{x} - \theta^{x}_0)    & (\theta^{y} - \theta^{y}_0) & (L - L_0) & (\phi - \phi_0)
                \end{bmatrix}$$
Where, $\mathbf{K}$, $\theta^{x}_0$, $\theta^{y}_0$, $L_0$ and $\phi_0$  is calculated from reference DNA while $\theta^{x}$, $\theta^{y}$, $L$ and $\phi$ is calculated for probe DNA from each frame.
This command gives output energy as a function of time in csv file and also average energies with error.
End of explanation
"""


%%bash

dnaMD globalEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 24 -estype ST -et "all" -gt "gmx analyze" -em "block"  -o energy_all_ST.csv

# Print first  and last 3 line of output file
echo "==========================================================="
echo "Deformation free energy of bound DNA as a function of time"
echo "==========================================================="
head -4 energy_all_ST.csv
printf ".\n.\n.\n"
tail -3 energy_all_ST.csv

"""
Explanation: Some of the plots from above data can be found here

Deformation free energy can be calculated as the following terms:

full : Use entire elastic matrix -- all motions with their coupling
diag : Use diagonal of elastic matrix -- all motions but no coupling
b1 : Only bending-1 motion
b2 : Only bending-2 motion
stretch : Only stretching motion
twist : Only Twisting motions
st_coupling : Only stretch-twist coupling motion
bs_coupling : Only Bending and stretching coupling
bt_coupling : Only Bending and Twisting coupling
bb_coupling : Only bending-1 and bending-2 coupling
bend : Both bending motions with their coupling
st : Stretching and twisting motions with their coupling
bs : Bending (b1, b2) and stretching motions with their coupling
bt : Bending (b1, b2) and twisting motions with their coupling

When all is used, all above terms were calculated.
End of explanation
"""


%%bash

dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 7 -o local_elasticity_4-7bps.csv

"""
Explanation: Same as above but only for stretching and twisting motions.
Local elastic properties
Local elastic properties can be caluclated using either local base-step parameters or local helical base-step parameters.
In case of base-step parameters: Shift ($Dx$), Slide ($Dy$), Rise ($Dz$), Tilt ($\tau$), Roll ($\rho$) and Twist ($\omega$), following elastic matrix is calculated.
$$
\mathbf{K}{base-step} = \begin{bmatrix}
K{Dx}        & K_{Dx,Dy}      & K_{Dx,Dz}      & K_{Dx,\tau}      & K_{Dx,\rho}      & K_{Dx,\omega} \
K_{Dx,Dy}     & K_{Dy}         & K_{Dy,Dz}      & K_{Dy,\tau}      & K_{Dy,\rho}      & K_{Dy,\omega} \
K_{Dx,Dz}     & K_{Dy,Dz}      & K_{Dz}         & K_{Dz,\tau}      & K_{Dz,\rho}      & K_{Dz,\omega} \
K_{Dx,\tau}   & K_{Dy,\tau}    & K_{Dz,\tau}    & K_{\tau}         & K_{\tau, \rho}   & K_{\tau,\omega} \
K_{Dx,\rho}   & K_{Dy,\rho}    & K_{Dz,\rho}    & K_{\tau, \rho}   & K_{\rho}         & K_{\rho,\omega} \
K_{Dx,\omega} & K_{Dy,\omega}  & K_{Dz,\omega}  & K_{\tau, \omega} & K_{\rho, \omega} & K_{\omega} \
\end{bmatrix}
$$
In case of helical-base-step parameters: x-displacement ($dx$), y-displacement ($dy$), h-rise ($h$), inclination ($\eta$), tip ($\theta$) and twist ($\Omega$), following elastic matrix is calculated.
$$
\mathbf{K}{helical-base-step} = \begin{bmatrix}
K{dx}        & K_{dx,dy}      & K_{dx,h}      & K_{dx,\eta}      & K_{dx,\theta}      & K_{dx,\Omega} \
K_{dx,dy}     & K_{dy}         & K_{dy,h}      & K_{dy,\eta}      & K_{dy,\theta}      & K_{dy,\Omega} \
K_{dx,h}      & K_{dy,h}       & K_{h}         & K_{h,\eta}       & K_{h,\theta}       & K_{h,\Omega} \
K_{dx,\eta}   & K_{dy,\eta}    & K_{h,\eta}    & K_{\eta}         & K_{\eta, \theta}   & K_{\eta,\Omega} \
K_{dx,\theta} & K_{dy,\theta}  & K_{h,\theta}  & K_{\eta, \theta} & K_{\theta}         & K_{\theta,\Omega} \
K_{dx,\Omega} & K_{dy,\Omega}  & K_{h,\Omega}  & K_{\eta, \Omega} & K_{\theta, \Omega} & K_{\Omega} \
\end{bmatrix}
$$
End of explanation
"""


%%bash

dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 7 -fgap 200 -ot local_elasticity_time_4-7bps.csv

# Print first 3 and last 3 line of first seven columns from output file
echo "====================================="
echo "Elastic matrix as a function of time"
echo "====================================="
head -4 local_elasticity_time_4-8bps.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, "... ... ..."}'
printf ".\n.\n.\n"
tail -3 local_elasticity_time_4-8bps.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, "... ... ..."}'

"""
Explanation: Local elastic properties as a function of time
Same command can be used to calculate elasticity as a function of time with option 
-ot/--output-time and save it in csv format file. This result can be used to check
their convergence. 

If this option is used, -fgap/--frame-gap is an essential option.
The output file is in csv format and can be opened as spreadsheet.

NOTE: Elastic properties cannot be calculated using a single frame because 
fluctuations are required. Therefore, here time means trajectory between zero 
time to given time.
End of explanation
"""


%%bash

dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -fgap 200 -gt "gmx analyze" -em "acf" -os local_elasticity_segments.csv

# Print first 3 and last 3 line of first seven columns from output file
echo "========================================"
echo "Elasticity as a function of DNA segments"
echo "========================================"
head -4 local_elasticity_segments.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, "... ... ..."}'
printf ".\n.\n.\n"
tail -3 local_elasticity_segments.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, "... ... ..."}'

"""
Explanation: Local elastic properties  of the consecutive overlapped DNA segments
Above method gives local elasticities of a small local segment of the DNA. However, we mostly interested
in large segment of the DNA. This large segment can be further divided into smaller local segments. 
For these smaller segments local elasticities can be calculated. Here these segments overlapped with each other.
Same command can be used to calculate elasticity as a function of time with option 
-os/--output-segments and save it in csv format file.

If this option is used, -fgap/--frame-gap is an essential option.
The output file is in csv format and can be opened as spreadsheet.
End of explanation
"""


%%bash

dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 10 -be 13 -et all -gt "gmx analyze" -em "block" -o local_energy_time_4-7bps.csv

# Print first 3 and last 3 line from output file
echo "==============================================="
echo "Local deformation energy as a function of time"
echo "==============================================="
head -4 local_energy_time_4-7bps.csv
printf ".\n.\n.\n"
tail -3 local_energy_time_4-7bps.csv

"""
Explanation: Local deformation energy of a local small segment
At first, elastic matrix from reference DNA (most often free or unbound DNA) is calculated 
and subsequently this matrix is used to calculate deformation free energy of probe DNA 
(most often bound DNA).
$$G = \frac{1}{2}\mathbf{xKx^T}$$
When helical='False'
$$\mathbf{K} = \mathbf{K}_{base-step}$$
$$\mathbf{x} =  \begin{bmatrix}
                  (Dx_{i}-Dx_0)  &  (Dy_i - Dy_0) & (Dz_i - Dz_0) & (\tau_i - \tau_0) &
                  (\rho_i - \rho_0) & (\omega_i - \omega_0)
              \end{bmatrix}$$
When helical='True'
$$\mathbf{K} = \mathbf{K}_{helical-base-step}$$
$$\mathbf{x} =  \begin{bmatrix}
                  (dx_{i}-dx_0)  &  (dy_i - dy_0) & (h_i - h_0) & (\eta_i - \eta_0) &
                  (\theta_i - \theta_0) & (\Omega_i - \Omega_0)
              \end{bmatrix}$$
End of explanation
"""


%%bash

dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 10 -be 13 -et all -helical -gt "gmx analyze" -em "block" -o local_helical_energy_time_4-7bps.csv

# Print first 3 and last 3 line from output file
echo "======================================================"
echo "Local helical deformation energy as a function of time"
echo "======================================================"
head -4 local_energy_time_4-7bps.csv
printf ".\n.\n.\n"
tail -3 local_energy_time_4-7bps.csv

"""
Explanation: Some of the plots from above data can be found 
here

Same as the above but energy is calculated using helical base-step parameters
End of explanation
"""


%%bash

dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -et all -gt "gmx analyze" -em "block" -os local_energy_segments.csv

# Print first 3 and last 3 line from output file
echo "=================================================="
echo "Local deformation energy as a function of segments"
echo "=================================================="
head -4 local_energy_segments.csv
printf ".\n.\n.\n"
tail -3 local_energy_segments.csv

"""
Explanation: Deformation energy of the consecutive overlapped DNA segments
Above method gives energy of a small local segment of the DNA. 
However, we mostly interested in large segment of the DNA. This large segment 
can be further divided into smaller local segments. For these smaller segments 
local deformation energy can be calculated. Here these segments overlapped with each other.
End of explanation
"""


%%bash

dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -et all -helical -gt "gmx analyze" -em "block" -os local_helical_energy_segments.csv

# Print first 3 and last 3 line from output file
echo "=========================================================="
echo "Local helical deformation energy as a function of segments"
echo "=========================================================="
head -4 local_helical_energy_segments.csv
printf ".\n.\n.\n"
tail -3 local_helical_energy_segments.csv

"""
Explanation: Some of the plots from above data can be found 
here

Same as the above but energy is calculated using helical base-step parameters
End of explanation
"""


