2018年11月14日 星期三
2018年11月12日 星期一
2018年11月11日 星期日
PDB FILE WRITER
Fortran:
! module, output a single frame (pdb "model") of the trajectory. This has | |
! nothing to do with OpenMM. | |
SUBROUTINE myWritePDBFrame(frameNum, timeInPs, energyInKcal) | |
use MyAtomInfo; implicit none | |
integer frameNum; real*8 timeInPs, energyInKcal | |
integer n | |
print "('MODEL',5X,I0)", frameNum | |
print "('REMARK 250 time=', F0.3, ' picoseconds; Energy=', F0.3, ' kcal/mole')", & | |
timeInPs, energyInKcal | |
do n = 1,NumAtoms | |
print "('ATOM ', I5, ' ', A4, ' SLT 1 ', 3F8.3, ' 1.00 0.00')", & | |
n, atoms(n)%pdb, atoms(n)%posInAng | |
end do | |
print "('ENDMDL')" | |
END SUBROUTINE |
static void | |
myWritePDBFrame(int frameNum, double timeInPs, double energyInKcal, | |
const MyAtomInfo atoms[]) | |
{ | |
// Write out in PDB format -- printf is so much more compact than formatted cout. | |
printf("MODEL %d\n", frameNum); | |
printf("REMARK 250 time=%.3f ps; energy=%.3f kcal/mole\n", | |
timeInPs, energyInKcal); | |
for (int n=0; *atoms[n].pdb; ++n) | |
printf("ATOM %5d %4s SLT 1 %8.3f%8.3f%8.3f 1.00 0.00\n", | |
n+1, atoms[n].pdb, | |
atoms[n].posInAng[0], atoms[n].posInAng[1], atoms[n].posInAng[2]); | |
printf("ENDMDL\n"); | |
}
|
參考
https://github.com/pandegroup/openmm/blob/master/examples/HelloSodiumChlorideInFortran.f90
https://github.com/pandegroup/openmm/blob/master/examples/HelloSodiumChloride.cpp
Lennard-Jones Fluid Properties
參考
2018年11月9日 星期五
Nose-Hoover NPT
Nose-Hoover chain length
現在分子模擬已經發展出許多thermostats和Barostats,最常用的莫非是利用Nose-Hoover演算法所建立的系統像是NVTNoseHoover和NPTMartynaTobiasKlein。這兩種方法都使用一種以上thermostats調控系統,每個thermostats彼此相互調節,使得溫度以及晶格震盪可以有效壓制,特別是模擬小系統效果非常好。Nose-Hoover chain length就是指串連幾個thermostats,chain length不影響平衡結果,但是平衡過程溫度震盪很大的話,可以增加chain length來減少溫度震盪。
參考
現在分子模擬已經發展出許多thermostats和Barostats,最常用的莫非是利用Nose-Hoover演算法所建立的系統像是NVTNoseHoover和NPTMartynaTobiasKlein。這兩種方法都使用一種以上thermostats調控系統,每個thermostats彼此相互調節,使得溫度以及晶格震盪可以有效壓制,特別是模擬小系統效果非常好。Nose-Hoover chain length就是指串連幾個thermostats,chain length不影響平衡結果,但是平衡過程溫度震盪很大的話,可以增加chain length來減少溫度震盪。
參考
2018年11月8日 星期四
tail corrections
在rc到無限遠的空間內散佈著i和j粒子,所以只是在算i和j粒子密度分佈。如果把這裡的密度思考成質量就會發現右項單位和utail單位不合。從講義上也可以確認rho的定義是"average number density"。
參考
http://www.iitg.ac.in/physics/fac/padmakumarp/Courses/PH442/PH442_Slides.pdf
根據Cameron F. Abrams的code會發現少一個N
但其實只是自己眼殘
參考
https://github.com/cameronabrams/codes/blob/master/src/mclj_npt.c
參考
http://www.iitg.ac.in/physics/fac/padmakumarp/Courses/PH442/PH442_Slides.pdf
根據Cameron F. Abrams的code會發現少一個N
但其實只是自己眼殘
參考
https://github.com/cameronabrams/codes/blob/master/src/mclj_npt.c
unit converter
Conventional units:
https://lammps.sandia.gov/doc/99/units.html
https://www.unitconverters.net/
Time
1 fs = 10^-15 s
1 ps = 10^-12 s
Pressure
1 pascal (Pa) = 1 newton/square meter (N/m2)
1 bar = 0.98692 atmosphere (atm) = 10^5 pascals (Pa)
1 atm = 101325 pascals (Pa)
1 pascals (Pa) = 0.0000098692 atm
1 GPa =10^9 Pa
1 MPa =10^6 Pa
Energy
1 joule (J) = 1 newton meter (Nm)
1 calorie (cal) = 4.184 joule (J)
R is related to the Boltzmann constant, k, by
R = k NA
where k = 1.3806 x 10^-23 J K^-1, and NA = 6.022 x 10^23 mol^-1
- distance = Angstroms
- time = femtoseconds
- mass = grams/mole
- temperature = degrees K
- pressure = atmospheres
- energy = Kcal/mole
- velocity = Angstroms/femtosecond
- force = grams/mole * Angstroms/femtosecond^2
- charge = +/- 1.0 is proton/electron
- distance = sigmas
- time = reduced LJ tau
- mass = ratio to unitless 1.0
- temperature = reduced LJ temp
- pressure = reduced LJ pressure
- energy = epsilons
- velocity = sigmas/tau
- force = reduced LJ force (sigmas/tau^2)
- charge = ratio to unitless 1.0
This listing of variables assumes conventional units; to convert to LJ reduced units, simply substitute the appropriate term from the list above. E.g. x is in sigmas in LJ units. Per-mole in any of the units simply means for 6.023 x 10^23 atoms.
參考Meaning Variable Units positions x Angstroms velocities v Angstroms / click (see below) forces f Kcal / (mole - Angstrom) masses mass gram / mole charges q electron units (-1 for an electron) (1 e.u. = 1.602 x 10^-19 coul) time --- clicks (1 click = 48.88821 fmsec) timestep dt clicks input timestep dt_in fmsec time convert dtfactor 48.88821 fmsec / click temperature t_current degrees K t_start t_stop input damping t_freq_in inverse fmsec internal temp t_freq inverse clicks damping dielec const dielectric 1.0 (unitless) Boltmann const boltz 0.001987191 Kcal / (mole - degree K) virial virial[xyz] Kcal/mole = r dot F pressure factor pfactor 68589.796 (convert internal to atmospheres) internal p_current Kcal / (mole - Angs^3) pressure p_start p_stop input press p_start_in atmospheres p_stop_in output press log file atmospheres input damping p_freq_in inverse time internal press p_freq inverse clicks damping pot eng e_potential Kcal/mole kin eng e_kinetic Kcal/mole eng convert efactor 332.0636 (Kcal - Ang) / (q^2 - mole) (convert Coulomb eng to Kcal/mole) LJ coeffs lja,ljb Kcal-Angs^(6,12)/mole bond various see force_fields file parameters 2,3,4-body terms
https://lammps.sandia.gov/doc/99/units.html
https://www.unitconverters.net/
Time
1 fs = 10^-15 s
1 ps = 10^-12 s
Pressure
1 pascal (Pa) = 1 newton/square meter (N/m2)
1 bar = 0.98692 atmosphere (atm) = 10^5 pascals (Pa)
1 atm = 101325 pascals (Pa)
1 pascals (Pa) = 0.0000098692 atm
1 GPa =10^9 Pa
1 MPa =10^6 Pa
1 joule (J) = 1 newton meter (Nm)
1 calorie (cal) = 4.184 joule (J)
Gas constant
pV = nRT R is related to the Boltzmann constant, k, by
R = k NA
where k = 1.3806 x 10^-23 J K^-1, and NA = 6.022 x 10^23 mol^-1
R with different units
8.31451 J K^-1 mol^-1
8.20578 x 10^-2 L atm K^-1 mol^-1
8.31451 x 10^-2 L bar K^-1 mol^-1
8.31451 Pa m^3 K^-1 mol^-1
1.98722 cal K^-1 mol^-1
1 Pa = 10^-3/4.184 kcal/m^3 = 10^-33/4.184 kcal/A^3 = 6.02*10^23/4.184/10^33 kcal/mol/A^3
1 GPa = 10^6 Pa = 6.02*10^29/4.184/10^33 kcal/mol/A^3 = 10^-4*6.02/4.184 kcal/mol/A^3
1 kcal/mol/A^3 = 4184*10^30/6/10^23 j/m^3 (Pa) * 0.0000098692 (atm/Pa) = 412927.328 atm
timestep (fs轉無因次)
t* = t * (epsilon / m / sigma^2)^1/2
t* = t([fs]) * [ sqrt( Kcal/mol * mol/g * Angstrom^-2 ) ] (1 Kcal = 4184J)
= t([fs]) * [ sqrt( 4184*Kg*meter^2/(second^2) * 1/g * 1/(10^-20 meter^2) ) ] (1 Angstrom = 10^-10 meter)
= t([fs]) * [ sqrt( 4184*1000 g*meter/(10^30 fs^2) * 1/g * 1/(10^-20 meter^2) ) ] (1 second = 10^15 femtosecond)
= t * 0.020455
pressure (atm轉無因次)
P* = P([atm]) sigma^3 / epsilon (1 atm = 101325 Pa = 101325 Kg/m/s^2)
= P([101325 Kg/meter/s^2]) * [ Angstrom^-3 * mol/Kcal ] (1 Kcal = 4184J)
= P([101325 Kg/meter/s^2]) * [10^-30 meter * mol*s^2/4184/Kg/meter] (1 mol = 6*10^23)
= P([101325 Kg/meter/s^2]) * [10^-30 meter * 6*10^23*s^2/4184/Kg/meter]
= P * 0.00001453
temperature (K轉無因次)
T* = T Kb / epsilon
= T([K]) * [ 1.38064852*10^-23 J/K * mol/Kcal ] (1 mol = 6*10^23)
= T([K]) * [ 1.38064852*10^-23 J/K * 6*10^23/4184 J ]
= T * 0.0019798975
1 Pa = 10^-3/4.184 kcal/m^3 = 10^-33/4.184 kcal/A^3 = 6.02*10^23/4.184/10^33 kcal/mol/A^3
1 GPa = 10^6 Pa = 6.02*10^29/4.184/10^33 kcal/mol/A^3 = 10^-4*6.02/4.184 kcal/mol/A^3
1 kcal/mol/A^3 = 4184*10^30/6/10^23 j/m^3 (Pa) * 0.0000098692 (atm/Pa) = 412927.328 atm
timestep (fs轉無因次)
t* = t * (epsilon / m / sigma^2)^1/2
t* = t([fs]) * [ sqrt( Kcal/mol * mol/g * Angstrom^-2 ) ] (1 Kcal = 4184J)
= t([fs]) * [ sqrt( 4184*Kg*meter^2/(second^2) * 1/g * 1/(10^-20 meter^2) ) ] (1 Angstrom = 10^-10 meter)
= t([fs]) * [ sqrt( 4184*1000 g*meter/(10^30 fs^2) * 1/g * 1/(10^-20 meter^2) ) ] (1 second = 10^15 femtosecond)
= t * 0.020455
pressure (atm轉無因次)
P* = P([atm]) sigma^3 / epsilon (1 atm = 101325 Pa = 101325 Kg/m/s^2)
= P([101325 Kg/meter/s^2]) * [ Angstrom^-3 * mol/Kcal ] (1 Kcal = 4184J)
= P([101325 Kg/meter/s^2]) * [10^-30 meter * mol*s^2/4184/Kg/meter] (1 mol = 6*10^23)
= P([101325 Kg/meter/s^2]) * [10^-30 meter * 6*10^23*s^2/4184/Kg/meter]
= P * 0.00001453
temperature (K轉無因次)
T* = T Kb / epsilon
= T([K]) * [ 1.38064852*10^-23 J/K * mol/Kcal ] (1 mol = 6*10^23)
= T([K]) * [ 1.38064852*10^-23 J/K * 6*10^23/4184 J ]
= T * 0.0019798975
2018年11月4日 星期日
GNU Scientific Library (gsl) on Mac
GNU Scientific Library是一種數值計算的C語言函式庫
https://www.gnu.org/software/gsl/doc/html/index.html
由於分子模擬的平衡速度為Maxwell-Boltzmann distribution,為了讓初始結構在平衡過程中不受到local minimum影響,可以利用gsl_ran_exponential設定系統初速度。
由於分子模擬的平衡速度為Maxwell-Boltzmann distribution,為了讓初始結構在平衡過程中不受到local minimum影響,可以利用gsl_ran_exponential設定系統初速度。
2018年11月3日 星期六
pybind11 for Mac
1. 三種方式安裝pybind11
$ brew install pybind11
$ pip3 install pybind11
$ conda install -c conda-forge pybind11
2. 寫example.cpp
參考
https://pybind11.readthedocs.io/en/stable/basics.html#header-and-namespace-conventions
3. 執行Mac的編譯指令產生example.cpython-36m-darwin.so
$ g++ -O3 -Wall -shared -std=c++11 -undefined dynamic_lookup `python3 -m pybind11 --includes` example.cpp -o example`python3-config --extension-suffix`
-undefined dynamic_lookup: Python gets linked dynamically without having all symbols resolved at link timepython3 -m pybind11 --includes: import pybind11 module
python3-config: output build options for python C/C++ extension
--extension-suffix: print the extension suffix used for binary extensions (產生example.cpython-36m-darwin.so)
參考
https://pybind11.readthedocs.io/en/stable/compiling.html#building-manually
https://osmocom.org/issues/1678
https://blog.csdn.net/IAlexanderI/article/details/81003225
https://helpmanual.io/man1/python3-config/
4. 利用python call c++ function
$ brew install pybind11
$ pip3 install pybind11
$ conda install -c conda-forge pybind11
2. 寫example.cpp
參考
https://pybind11.readthedocs.io/en/stable/basics.html#header-and-namespace-conventions
3. 執行Mac的編譯指令產生example.cpython-36m-darwin.so
$ g++ -O3 -Wall -shared -std=c++11 -undefined dynamic_lookup `python3 -m pybind11 --includes` example.cpp -o example`python3-config --extension-suffix`
-undefined dynamic_lookup: Python gets linked dynamically without having all symbols resolved at link timepython3 -m pybind11 --includes: import pybind11 module
python3-config: output build options for python C/C++ extension
--extension-suffix: print the extension suffix used for binary extensions (產生example.cpython-36m-darwin.so)
參考
https://pybind11.readthedocs.io/en/stable/compiling.html#building-manually
https://osmocom.org/issues/1678
https://blog.csdn.net/IAlexanderI/article/details/81003225
https://helpmanual.io/man1/python3-config/
4. 利用python call c++ function
訂閱:
文章 (Atom)