Differences

This shows you the differences between two versions of the page.

--- surface_area_of_lipid [2021/05/25 13:36] – [Overall Average SA/Lipid] admin
+++ surface_area_of_lipid [2026/09/20 00:36] (current) – refactored_scripts edit
@@ Line 3: / Line 3: @@
 ==== Overall Average SA/Lipid ====
-(authors: Xiaohong Zhuang & Jeffery Klauda)
+(authors: Xiaohong Zhuang & Jeffery Klauda)\\
+Refactored by Eric Zhang (2026)
 === Introduction ===
-This simple analysis is used to estimate the overall average surface area (SA) per lipid. If you are analyzing a single component membrane then this is equivalent to the component's SA/lipid. Otherwise, this is the average of all lipids.
+This simple analysis is used to estimate the overall average surface area (SA) per lipid. If you are analyzing a single component membrane then this is equivalent to the component's SA/lipid. Otherwise, this is the average of all lipids. UPDATE: We highly recommend using the modern Python scripts for analyses (instructions listed at the very bottom). Also, those using the modern scripts DO NOT need to run the Equilibrium Test and Area Compressibility Analyses separately (this is automatically handled).
-Those at UMD should use the DT2 path: **/lustre/jbklauda/scripts/area_1_lipid**
+Old scripts:\\
+Those at UMD should use the ZT1 path: **/afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/area_1_lipid/**\\
+Outside of UMD a sample script is given in the gzip file: {{ :area_1_lipid.gz |}}\\
+Example of analysis result in an Excel File: {{ :e.g._sa_1_lipid.xlsx |}}\\
-Outside of UMD a sample script is given in the gzip file: {{ :area_1_lipid.gz |}}
+Modern Python Scripts:\\
+{{ ::area_1_lipid_v2.tar.gz |}}
-Example of analysis result in an Excel File: {{ :e.g._sa_1_lipid.xlsx |}}
+=== Old Script Instructions ===
-=== Instructions ===
 . Copy the folder //area_1_lipid// from the path highlighted above to your directory has contains all dyn files. //area_1_lipid// should become a subdirectory of your folder with the dyn*.dcd files.
@@ Line 36: / Line 39: @@
   * area2.scr is used to obtain the areas from the the first dyn file (likely dyn1.xst) to the end
-  * The //awk// line creates of time series of the area assuming time is in fs
+  * The ''awk'' line creates of time series of the area assuming time is in fs
-  * the //sed// line reduces the printing of data (every 50 lines) in a reduced file...this may not be needed if the printout is not frequent.
+  * the ''sed'' line reduces the printing of data (every 50 lines) in a reduced file...this may not be needed if the printout is not frequent.
-You will need to run the scr file with type the following on the command line: //./1_run_area.scr//
+You will need to run the scr file with type the following on the command line: ''./1_run_area.scr''
 After run **1_run_area.scr**, the files //area_time_acc.dat// and //area_time_acc_excel.dat// will be generated. For //area_time_acc.dat//, the first column is time in nanosecond (ns), the 2nd column is the surface area per lipid at each time (step), and the 3rd column is the cumulative average of surface area per lipid.
-.2. Generate the plot of area vs. time in order to find the dyn that the system reaches equilibrium
+.2. Generate the plot of area vs. time in order to estimate the time the system reaches equilibrium with respect to the SA/lipid.
-.2.1. If you are using WinSCP, you may run the script by command:  ''2_plot_sa.scr'' to use Gnuplot to generate plot //2_sa.tif//. And then you can open and view the image in WinSCP by right-click on the file. Plot using Gnuplot I personally think is convenient in the linux system but it is optional. An example of the plot is shown below.
+.2.1. Using GnuPlot on LINUX:  ''./2_plot_sa.scr'' command is to use Gnuplot to generate plot //2_sa.tif//. An example of the plot is shown below.
 {{:fig1exsaplip.png?200|}}
 Figure 1. Example surface area per lipid plot
-.2.2. If you don’t use WinSCP,or if you are not comfortable with Gnuplot, you can also use Excel to plot //area_time_acc_excel.dat// in Excel. The example of the plot in Excel is attached.
+It is recommended to adjust the GnuPlot script to adjust the axis scale to best represent your data. Using nano or vim as an editor of //2_sa.gnu//
-Based on the area plot, we can determine when (at which dyn) the system reaches equilibrium.
+<code>
-The equilibrium dyn number is the dyn file number that your cumulative (or running) average surface area per lipid (sa/lip) (red data points in Fig.1) becomes constant for around 20-30ns (i.e 10-15 dyn, refer to note). Since the running average including the beginning high values,we will also determine equilibrium based on the area per lipid (blue data points in Fig.1) that fluctuate around the same center line).
+# show margin
+plot [0:100] [35:50] \
+'area_time_acc.dat' using 2:3 with point  pt 2 ps 2 lw 4 lc rgb 'blue' title '{/DejaVuSerif-Bold=45 Area/Lip.}', \
+'area_time_acc.dat' using 2:4 with point  pt 6 ps 3 lw 4 lc rgb 'red'  title '{/DejaVuSerif-Bold=45 Cm. Avg.}'
+</code>
-e.g. For in Fig.1, the cumulative average area (red) indicate that system reach equlibrium at around 30ns, and the areas (blue) fluctuate around the same center line, so we use the data from 50-100ns (dyn 26-50) or 60-100ns(dyn 31-50)  to calculate the average area/lip.
+The ranges after the //plot// command can be adjusted for the time (in ns) and SA. The other lines can be adjusted to format personally desired.
+.2.2. You can also use Excel to plot //area_time_acc_excel.dat// in Excel. The example of the plot in Excel is attached.
-.  Based on equilibrated dyn range on step 2, select the data points from equilibrated dyn to the end to calculate the block-average surface areas.
+.3. Estimating time for equilibration.
-.1 In 3_cal_avg_std.scr, please update the total number of last dyn*.dcd files to calculated the averages.
- <code> set Ndyn = 20 </code>
-The example states that you want to do the averages from the last 20 dyn*dcd files. Let's say you have 50 dyn*.dcd files, then the first 30 files are ignored. If each file had 2ns of data, then you are calculating the average from 60-100ns.
+Based on the area plot, you can determine when (at which dyn) the system reaches equilibrium.
+The equilibrium dyn number is the dyn file number that your cumulative (or running) average surface area per lipid (sa/lip) (red data points in Fig.1) becomes constant (i.e 10-15 dyn, refer to note). Since the running average including the beginning high values, we will also determine equilibrium based on the area per lipid (blue data points in Fig.1) that fluctuate around the same center line).
+For in Fig.1, the cumulative average area (red) indicate that system reach equlibrium at around 30ns, and the areas (blue) fluctuate around the same center line, so we use the data from 50-100ns (dyn 26-50) or 60-100ns(dyn 31-50) to calculate the average area/lip. You will need to do the same for your system.
-# This following command calculates the block average, the output file name (//area_avg_std.dat//) and the block size (1000) in time steps i.e. each dyn is treat as a block
+. Obtaining Statistics of SA/Lipid
-.2 Calculate the block-average surface areas by command: ''3_cal_avg_std.scr''
+Based on equilibrated dyn range determined in Step 2.3, select the data points from equilibrated dyn to the end to calculate the block-average surface areas.
-# Store the output file name (//area_avg_std.dat//) and the block size (# of steps in a DCD file=Ndata) in time steps into a file //area-block.inp//, which is the input data file for //area-block.exe//
+.1 ''3_cal_avg_std.scr''
-<code> echo area_avg_std.dat $Ndata > area-block.inp
-area-block.exe < area-block.inp </code>
-Note: For  //#define MAX_R 100000//  in area-block.c, you might have to change the number to a larger one (larger than your # data points) and recompile it.
-As you can see, the output average area data are saved in area_avg_std.dat. The result is given in average+- population standard deviation.
+Update the total number of last dyn*.dcd files to calculated the averages.
-. Calculate averages, sample standard deviation, and standard error of block averages by command: ''4_run_calc_stderr.scr''
+ <code> set Ndyn = 20 </code>
-The file //avg_stderr.dat// will be generated that will show the final values that we need, which are the average of the block averages and the standard error of the block averages, i.e. SA_avg/lip +- Standard error. You may refer to the attached Excel file which also demonstrate how the standard error is calculated if you are interested.
+The example states that you want to do the averages from the last 20 dyn*dcd files. Let's say you have 50 dyn*.dcd files, then the first 30 files are ignored. If each file had 2ns of data, then you are calculating the average from 60-100ns.
-Note:
+.2 Calculate the block-average surface areas by command: ''./3_calc_avg_std.scr''
-The corresponding nanosecond are based on the "number steps" and "timestep" that we set.
-For example, in the //dyn-2.inp// file, we see  there a code line:  "numsteps            1000000           ;
-# run stops when this step is reached"
-and another code line:  <code> timestep            2.0  ;# fs/step </code>
-So you multiple 2fs/step  with 1,000,000  steps per dyn file,  which makes is 2ns per dyn file. (since 1fs=10^(-15) s, 1ns==10^(-9) s )
-==== Get the equilibrium data (test) ====
+**Note: **For  //#define MAX_R 1000000//  in area-block.c, you might have to change the number to a larger one (larger than your # data points) and recompile it.
-(author: Yalun Yu,  email: alanyu17@terpmail.umd.edu / yalun.research@gmail.com)
+As you can see, the output average area data are saved in //area_avg_std.dat//. The result is given in average+- population standard deviation.
-< Email me if anything is wrong >
-This is a brief instruction to obtain the time range of equilibrium (and also the minimum block size can be used for standard error calculation). Based on these, calculation of SA/lipid and Ka is done automatically for your lipid only system. Run Xiao's //1_run_area.scr// before doing this. Please put the scripts in the same folder as you run //1_run_area.scr//. You will need pymbar package (in python) for this:
+.3 Standard Errors from Block Average
-. Install Anaconda (https://www.anaconda.com/distribution/).
+Calculate averages, sample standard deviation, and standard error of block averages by command: ''./4_run_calc_sterr.scr''
-   a) Download the python3.7 version for Linux [I used the '64-Bit (x86) Installer (652.5 MB)']
+The file //avg_stderr.dat// will be generated that will show the final values that we need, which are the average of the block averages and the standard error of the block averages, i.e. SA_avg/lip +- Standard error. You may refer to the attached Excel file which also demonstrate how the standard error is calculated if you are interested.
-   and put the .sh installer in your home directory (or any other place you can access and execute it) on DT2.
-   b) make the .sh file an executable (chmod u+x $NameOfInstaller.sh)
-   c) Install by running "./$NameOfInstaller.sh", keep in mind where you have installed it.
-. Make sure the Anaconda installation path is included (see below) in your ~/.cshrc.mine (if using csh) or ~/.bashrc.mine (if using bash) and source the .mine file.
-   In ~/.cshrc.mine, add this line:
-   set path = ($HOME/anaconda3/bin $path)
-   or in ~/.bashrc.mine, add this line:
-   export PATH="$HOME/anaconda3/bin:$PATH" [may differ depending on where Anaconda is installed, see 1(c)]
-. create a conda environment and activate the env. (see conda's website for instructions)
-. Install the pymbar package.
-   conda install -c omnia pymbar
-Download the script: {{ :area_ka.tar.gz |}}
-Change the temperature and number of lipids per leaflet in //pymbar_area_ka.scr//:
+==== Modern Python Script Instructions ====
-   set temp=323.15
+These scripts simplify and automate many of the steps in the old scripts. They also run ~10-50x faster. As mentioned earlier, they automatically check equilibrium and derive area compressibility, so you don't need the separate analysis listed elsewhere in the wiki. If there are any questions, please reach out to ezhang32@umd.edu.
-   set nlip=36
-Change the initial guess (or your preference) for block size (in ns) and the steps size (in ns) of your output (usually 1ps, but make sure to match your xst freq):
+=== Step 1: Edit ''config.py'' ===
-  area_handling.py 10 0.001 >& all.dat
+Open ''config.py''. Read the comments, then specify the number of lipids per leaflet and the path to your dyn/namd folder with dyn/dcd files.
-After running //pymbar_area_ka.scr//, you'll see //all.dat//, SA/lipid is in angstrom^2 and Ka is in N/m.
+Note: If you have an asymmetric bilayer, just let nlip = # of lipids in the top leaflet, finish the rest of the instructions, then change nlip = # of lipids in the bottom leaflet, then rerun the scripts. You'll realize that Equilibrium and Ka results don't depend on nlip.
+=== Step 2: Run Scripts ===
+Once config.py is populated, run ''./0_setup.py''. If you see error messages, just go back and fix config.py as directed and try rerunning. If system_info.txt is successfully generated, run ''./1_run_area.scr''. Once this finishes, a new results directory will appear, which contains plots and equilibirum / ka results.