網頁

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

C++:
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

參考

The purpose of these pages is to provide some explicit results from molecular dynamics and Monte Carlo simulations for the Lennard-Jones fluid. It is intended to provide guides for testing codes. Reproducing these results is a test of the correctness of codes, either written by the user or obtained elsewhere. The explicit conditions for each of the sets of results are supplied so that meaningful comparisons of your results with the ones listed here are possible.

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來減少溫度震盪。

參考

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

unit converter

Conventional units:
  • 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
LJ reduced units:
  • 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

Energy
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

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設定系統初速度。

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

















$ python3 example.py