To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 639ee84a authored by dkammer's avatar dkammer
Browse files

week 7: group 2 added

parent d9e89b02
build
\ No newline at end of file
cmake_minimum_required(VERSION 3.1)
project(tehpc)
enable_language(CXX)
# flags for compilation (options)
set(CMAKE_CXX_FLAGS "-Wall -Wextra -pedantic -std=c++11 " CACHE STRING "" FORCE)
set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG " CACHE STRING "" FORCE )
# alternative to add_library (but with option)
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.")
mark_as_advanced(BUILD_SHARED_LIBS)
# build type
set(CMAKE_BUILD_TYPE "${CMAKE_BUILD_TYPE}" CACHE STRING
"Build options: None Debug Release RelWithDebInfo MinSizeRel."
FORCE )
# include directories
set(TEHPC_INCLUDE_DIRS
"${CMAKE_CURRENT_BINARY_DIR}/src"
"${CMAKE_CURRENT_SOURCE_DIR}/src"
)
include_directories(${TEHPC_INCLUDE_DIRS})
find_package(Armadillo REQUIRED)
include_directories(${ARMADILLO_INCLUDE_DIRS})
# source directory
add_subdirectory(src)
# ------------------------------------------------------
# examples
# define add_simulation function
function(add_simulation SIM_EXE)
add_executable(${SIM_EXE} ${SIM_EXE}.cpp)
target_link_libraries(${SIM_EXE} tehpc)
endfunction()
option(TEHPC_EXAMPLES "examples" OFF)
if(TEHPC_EXAMPLES)
add_subdirectory(examples)
endif()
This diff is collapsed.
add_simulation(basic_example)
configure_file(basic_example_1.inp basic_example_1.inp COPYONLY)
configure_file(basic_example_1.rest basic_example_1.rest COPYONLY)
configure_file(basic_example_1.run basic_example_1.run COPYONLY)
#include <iostream>
#include "Model.hpp"
int main(int argc, char *argv[])
{
std::string fname;
std::string odir;
if(argc<3) {
std::cerr << "Not enough arguments:"
<< " ./basic_example <sim-name> <output-directory>" << std::endl;
return 0;
}
else {
fname = argv[1];
odir = argv[2];
}
Model my_model;
std::cout << "read input" << std::endl;
my_model.read_input(fname+".inp");
std::cout << "read restart" << std::endl;
my_model.read_restart(fname+".rest");
std::cout << "init" << std::endl;
my_model.init();
my_model.init_dump(odir, fname);
my_model.dump();
std::cout << "start time stepping" << std::endl;
int nb_increments = my_model.get_nb_time_steps() - my_model.get_time_step();
for(int i=0; i<nb_increments; ++i) {
my_model.increment_time_step();
my_model.velverlet();
my_model.compute_ekin();
my_model.dump();
}
std::cout << "end time stepping" << std::endl;
return 0;
}
108 # natoms
39.948 # mass in AMU
0.2379 # epsilon in kcal/mol
3.405 # sigma in angstrom
8.5 # rcut in angstrom
17.1580 # box length (in angstrom)
10000 # nr MD steps
5.0 # MD time step (in fs)
100 # output print frequency
6.67103294321331 -10.6146871435653 12.6336939877734
1.06574058650169 -3.33432278188177 -2.59038677851747
-1.78412295775301 -16.5259458407765 4.61680014503288
9.81030288841332 -5.82405963687712 3.65645974729484
15.3247211142637 -4.78053306974291 8.26672732196501
2.11253838645732 -4.72955953073222 6.75738627120688
-4.56209227856753 -3.20329426260617 5.55255473461665
5.86663918205102 -10.4139411857802 1.86416360433577
5.21527405147765 -28.9277726235844 6.19270143781154
13.7051046236876 -4.97120493339350 1.94629699552930
4.19613239880890 -14.9790822892036 16.5333880558744
16.7521502712298 -10.0804241645332 14.5028498925954
12.6350229329042 -1.40635636262803 2.38627163152721
-1.08439966554446 2.41557722385393 8.37791671838088
8.96146834305310 -19.0272834036820 4.05806985477701
11.1459036450055 3.18552372924961 6.723294858788574E-002
16.0690989728094 -1.03744620980343 18.5655733607861
18.6984479366938 -17.6491643729990 3.60417710784249
-3.13161415889863 -6.75371567663538 5.00611042020118
17.4072659061416 -7.88472086048433 7.56603710600785
19.6499870844957 -10.2813545329547 3.47269993677700
9.14358801056816 10.5769423320094 11.0434815971320
12.3852978120607 -11.0158733347720 20.0337606688683
6.18023614761844 -4.34746306558901 5.42729289210231
0.436403345315857 -1.36655907890507 7.37333479343561
-0.920905510346089 -3.30961460032849 21.5616988921137
3.31118739595332 -20.7760460532019 20.3172712920931
0.507105497066932 -21.6567949737129 0.933840061914709
14.4213709093942 -20.3431287076219 15.4309530884964
5.56191131389967 -7.25872501357231 12.4568959982054
3.06146523952825 -14.4087233421704 6.69325558343884
-1.37943237515652 -10.1053635146924 4.64502788848255
7.52415338992668 -14.4733254864890 18.7673317352093
-5.87703921849131 1.37272405947699 13.4338857916159
0.830260163129983 5.33087087306913 11.0816141860704
2.70133951013669 -18.7636102294181 0.120055694867087
9.27050983701818 -10.9072106043863 0.889331660756358
9.08594738166258 -17.5221802255344 0.686386365370126
5.25139950502609 -17.4969690646468 19.3171785990622
13.2681025537439 -13.8523638332393 5.66497231845133
15.2101790729727 2.32029736635329 11.6969471206408
11.2744379312723 -3.31524720571622 9.80174695909147
11.3672221692278 -7.00573339367229 -8.863893204687781E-002
12.6807143612004 -0.189603658778157 16.3325462080353
13.1935341627820 -1.62543720726316 29.5955788172380
15.3745579360649 -22.0045244935217 -4.92637554192936
16.4280280206508 -8.26617967392137 11.4518786798574
17.8237371703412 -15.2862663598538 17.6996045132077
13.0641438321788 -10.8486456984713 16.0847503564452
7.73638629377528 1.84526543149574 15.1833781698611
10.2571128578188 -11.9445746958196 14.2192606487898
3.77699860251963 2.92447970812240 2.90610322685796
9.87198946567047 -1.85194776144305 14.3336221929169
9.34121495781715 -0.110266095177049 10.2421994884800
6.42953094552702 -11.9658267479721 -1.27607222782124
0.312075928567517 -13.5209321587520 4.30576572830279
21.8849306695113 -8.32941966975701 15.7831238633949
-1.73962014751743 -6.92520932653883 15.0877890024375
15.3258711078380 -10.9370636474509 8.75014860388744
5.62898476019133 1.22206410978867 9.13368454732895
4.36630448913605 -7.17908827942928 2.46370171331934
6.36905055137017 2.140177685278017E-002 12.4806258246867
-8.24590606548181 -8.40492557103603 14.5774496555802
6.39507074174614 -8.93453582347402 9.13184152048137
3.17311027831434 -10.4109670230052 13.4136189406149
8.55871765808877 -11.5146971263885 4.70680651719831
3.77957617856256 -5.11519004071282 -0.633083143544725
3.82062013524474 -1.20262776483639 6.23531431086906
4.25706748939141 -3.63430271674530 13.2690508392963
1.92348687154464 5.96702632661158 7.62724137297314
12.0051777788663 -13.8263585530837 9.18079560803348
15.0453466485826 -14.2605922098574 15.3530811864134
24.4805821494876 -4.72207741599974 -1.99596985109394
19.4767646609806 1.77442931216133 10.3617496397394
11.9756015797765 -5.07615768419422 13.6907379724847
13.9394476435385 5.54214528144574 -5.00708130293960
17.2465568292747 -1.41716902620044 10.9507552983048
16.5494040530701 -0.264088515793768 15.1391134707549
8.16231980799074 9.84631749751864 0.796210036557643
8.40323357518367 3.35027973039626 11.5913687598978
23.3498272041744 -1.40692495029280 15.9949128751187
7.23312060426592 -2.05938049655674 7.98702327070333
1.04256284786092 -8.15292748227513 17.2074386421529
3.01310260723471 -11.5097160369980 16.8994512032187
6.88649054507511 -8.15930685942754 4.82906099336082
5.06237451942686 -12.2934797858039 9.79914937310931
16.6129272937896 -11.8133622739936 18.3417827366773
2.69214517473574 -8.28564119894275 10.4601950916785
2.79781292393336 0.114875238778070 13.8206239519260
7.62268397733944 -3.17673280498871 -5.53083875591493
1.50744766634734 -13.8205771069070 13.9164676875935
6.67307640084798 -16.0925383485583 5.25789452573394
10.4390555563715 -16.6707469563248 6.91799417023151
5.62891628767610 -5.22859253372635 9.43277399195221
14.1634497645406 -14.5361596069532 19.0790815297552
10.6505276709263 -3.69963830907047 17.6508669507771
3.70250355408209 -7.94643191534542 6.34294974462814
8.29060739488461 3.39283070761242 24.9724714610759
6.79682924336071 -20.9921316604048 1.41649721746698
19.1829730396108 -6.67133198675445 13.8686185773761
9.15960737692621 11.8185252568240 24.7562364347319
22.0227224042279 -13.7415388017328 13.2127365927663
18.6458829445424 -4.67218013657984 10.5730480824629
4.08992423327782 -2.00850424827988 9.63570479533125
14.6660200866878 -8.57920230381516 18.6426266424536
13.7138121972980 -0.358709470051391 8.46443111624972
18.2065606702052 -7.01462311401517 21.0114091185902
-6.65121935702062 1.49356278110318 20.6354475761263
-1.5643224621482283e-03 4.8497508563925346e-04 -4.3352481732883966e-04
4.1676710257651452e-04 2.2858522230176587e-05 -6.1985040462745732e-04
-7.5611349562333923e-04 4.0710138209103827e-04 -4.6520198934056357e-04
-1.0716463354369357e-03 1.9399472344407333e-03 3.6805207892014988e-05
1.8845309488827741e-03 -9.6929654979601939e-04 -5.0590631889323786e-04
6.3850936177063688e-04 2.9439961625943489e-03 2.0290627741083355e-04
1.5417020528064903e-05 1.7854479210267119e-05 1.0061778727049613e-03
-1.4096631765204916e-03 6.6920838746400600e-04 -1.1440332398881235e-03
-1.1623287197336312e-03 1.6266211340761735e-03 1.3040008893057014e-03
1.3961661508013856e-03 3.5480038145712736e-05 -9.9101224232255817e-04
-7.8119282252798756e-04 8.0770620385540233e-04 -1.4741616364769391e-04
-8.4465150294226528e-04 -1.9923244287298695e-03 -6.2089699425759077e-04
4.8199317804818579e-04 2.6851883010136799e-04 2.8946202010011360e-05
-8.2433856370986945e-04 2.8766236402950288e-04 -6.7726146942851230e-04
1.0759225321624727e-03 -5.0652068703702471e-04 -8.1558579057040967e-04
-4.1654885401277715e-04 -2.1963115969871258e-03 4.7173893044170949e-04
3.3822271329348960e-03 -1.1094925187093363e-04 1.7960828583444049e-03
-4.5419241078412093e-04 -5.9345998636867009e-05 -5.7258946451372564e-04
-1.1899036421980488e-03 9.3983719885128869e-05 -5.8462360685180778e-04
3.7172749436126885e-04 1.9340413720264325e-03 1.4034788302775946e-03
1.7974246348577072e-04 1.0193612025409288e-03 4.9370272372543282e-04
-9.4461498382508897e-04 -1.5026896525989781e-03 9.1479102888999955e-05
-1.8939664944843942e-03 5.6369115520127778e-04 2.4325248949145196e-04
5.8057892769663175e-04 1.2299831371874594e-04 2.9629659010955246e-04
8.3145733665166094e-04 1.6024783547625206e-03 2.2561613876544880e-03
-3.3459474408059099e-04 -9.0138124644271826e-04 4.6451973781908732e-05
-1.4661492298318566e-03 -9.6700324227027277e-04 -1.2402619350327583e-03
2.2294716528019855e-04 -7.1757354452299765e-04 -1.3359945828925550e-03
3.5744852428510476e-04 -4.7176712674574758e-04 -1.2299995762011908e-03
-7.1013397780573531e-04 2.1818976630239278e-03 3.0517895276767340e-03
7.9633427864940574e-04 -3.3213284823899663e-05 -8.8375819004122126e-04
-1.0048962733085658e-03 -1.0338986950231507e-03 -1.2128201825351069e-03
6.1036336377894929e-04 9.2701621970857817e-04 -1.2537972128693620e-03
1.0520429442651071e-03 -8.4866067556637161e-04 -1.3560525458749155e-05
2.5814213260471724e-04 -7.5314897639645743e-04 8.0713082892823452e-04
-3.6385410005198736e-03 8.8146280046642606e-04 5.6477118870225215e-05
9.6333620309687013e-04 1.7133724036878372e-04 -6.2326101449483172e-04
1.3552473384443623e-03 6.8479283464214743e-04 1.7273048455913563e-03
2.5884589579386447e-03 -1.4321360648116401e-03 -3.0829116871650114e-04
-1.9464964022266687e-03 -6.1160752453349304e-04 1.2508357692460391e-03
2.0354460092287877e-03 -8.0047183425335712e-04 -7.6779326062273879e-05
-8.8078983167544410e-04 2.0527738907615478e-03 -4.8490967626695514e-04
-1.6226833552402415e-03 -4.4509205363786440e-04 5.7985161430016720e-04
-1.3846623789497004e-03 -8.1591443760356856e-04 3.5239524960511482e-04
8.4581719582112696e-05 4.4027379945299169e-04 7.8825448363930943e-04
1.9101264956063397e-04 -3.7541323533216431e-05 2.0232815584298586e-04
1.0507771398316642e-03 2.1703338837872421e-06 1.8721233841267814e-03
-1.1764360038119695e-03 -2.7051268798865438e-03 4.2626857615829489e-04
1.4811895380049227e-03 -8.6900359536064796e-04 -1.3508346585503553e-03
1.3876291076744002e-03 6.5070110767351664e-04 -6.1537705602631544e-05
-1.5265028980997747e-04 3.5305182092379836e-03 -3.5107439080411313e-04
1.2593135527188843e-03 -3.4770110188951111e-04 -3.2219440087892302e-03
-4.8839641798958239e-04 -1.0781245860584732e-03 -7.6447169214134084e-04
2.9537856402441307e-03 3.5708460303038737e-04 4.3965885540392521e-04
-9.2354805883572888e-04 -1.3750879594929264e-04 1.7566816252148881e-03
-8.4649245109681394e-05 -2.5830032329337755e-03 6.2136732242066585e-04
-7.8863228626065347e-04 -4.7987306164718428e-05 -7.6918205615324461e-05
-7.5715366119523408e-04 -8.8223969491051666e-04 -1.1130209071318147e-03
-1.5829995718074550e-03 -1.0023130948058853e-03 1.0105693656524083e-03
-9.7237564644641694e-04 -4.2306183690077959e-03 -1.7993911047256889e-03
2.5009355760304557e-03 -6.9762487866167910e-04 1.4240088170272963e-03
-2.1973523929318302e-04 2.2092313541290937e-05 8.3472986487536491e-05
1.1022471068800079e-03 -3.1173869715215023e-03 2.2862705719842859e-03
-1.2765414792041748e-03 8.2954964873239924e-04 8.6393605486046691e-07
1.8024567462803108e-03 -3.2357868878445974e-04 7.2049328401194016e-04
1.3369574268302883e-03 2.4197150221213766e-04 -2.4078296952019445e-04
-6.4367273480117029e-04 9.4544158572617176e-05 1.2224424402153017e-03
7.3861349744237207e-04 6.1538146539837915e-04 -1.0659904512438458e-04
-1.3869668138548187e-03 -1.5135537737600894e-03 -1.7700632516888342e-03
-1.2783209148195544e-03 -6.4363699635212640e-04 -6.3784434948207917e-04
-1.1486989458051769e-03 2.2414503838807933e-03 5.7721354683718497e-05
1.3168895876064066e-03 1.9570675664972955e-03 4.0936567386207711e-04
6.3073775671656379e-04 6.9398368370810696e-04 -2.7072142801711462e-03
-3.2393828021454760e-04 1.8032188549262072e-03 -8.7265713577056329e-04
2.0422470188389682e-04 -7.7369026089478630e-04 -1.3328071442246694e-03
-2.2226907938257076e-05 4.7885157053263603e-04 6.0068136092359822e-05
-1.4298413641215028e-03 -1.1674283659453153e-04 -3.7073258095465290e-04
-2.0611289910871505e-03 -7.1332431440348544e-04 1.2809266282339321e-04
-7.4275381512506605e-04 1.8409280327853588e-03 -1.4613152937356479e-03
-1.1017089795515253e-03 9.2440680189215501e-04 4.1040985295765773e-05
5.3072557614502406e-04 -1.9182563618355668e-03 6.6929180130627429e-04
-1.6469389079891580e-03 -3.3913862865732141e-04 -3.4603537121251295e-04
-7.4215102626457323e-04 -9.8513589889560684e-05 8.8892011842225927e-04
-3.0986258142986708e-04 2.3322786412505996e-03 1.1214740271337398e-03
-2.2588473708594786e-03 -3.7370321664764982e-04 2.6199219861683513e-03
1.5543286843928208e-03 3.0860432881096969e-04 -2.3276316538484460e-03
1.4014253099048656e-03 1.6989378445509686e-04 1.0862067754042389e-04
1.5485169419261718e-03 -2.5482574014883033e-04 1.6549301081433053e-03
8.3679284925556360e-04 2.6544605630854541e-03 -3.1978397990891830e-04
2.8867805091229137e-04 -9.5105410949532254e-04 -4.1906535931744601e-04
1.5677334410636003e-03 3.8398276164589067e-04 5.4230186211901119e-04
1.8828375897567346e-03 2.5974363433841984e-04 -5.4853918242744198e-04
-7.2191331833469923e-04 -9.3121561995775873e-04 -1.0102462549550035e-03
-2.5645967264587824e-04 -3.8300695089345021e-04 7.5308084406414243e-06
-6.9678200159030467e-04 1.9053736728726351e-03 -5.5375918909208023e-04
1.8765197796989418e-03 2.2251752214693724e-04 -8.5240569555193027e-04
6.5009741785571464e-04 1.5502364257396865e-03 -1.5320902037263414e-03
1.1552195131152949e-03 5.2576645590408585e-04 -9.7934306024165957e-05
-1.4419828777329281e-03 -2.0581943028127967e-03 -6.1455516747372426e-04
-1.8081813454125322e-03 2.0392408575730525e-04 5.2891537053608253e-04
8.3338854452045131e-04 8.4885016567212471e-05 -1.0314730467739900e-03
1.7861843205888382e-03 -1.4959403377944355e-03 5.0308057595438297e-04
-6.0451783799718114e-04 -1.9369920772084048e-03 8.3982235283629170e-04
8.7324422742439877e-04 1.6025909516322815e-04 -6.4317473328755128e-04
-4.5793392919712261e-04 -2.1744186189341037e-04 2.5912834893712798e-03
1.5850650475339276e-03 -6.4434145921806031e-06 -5.7440449281686960e-04
-2.0204313429923672e-03 6.3428280954377821e-05 2.5115932933766497e-03
-7.6107091738481837e-05 6.5742688995404971e-04 -2.9040212591161766e-04
odir=~/dev/output
name=basic_example_1
./basic_example $name $odir 2>&1 | tee $name.progress
cp $name.inp $odir/$name.inp
cp $name.run $odir/$name.run
cp $name.progress $odir/$name.progress
set(TEHPC_HEADER
Model.hpp
#Mesh.hpp
)
set(TEHPC_SRC
Model.cpp
#Mesh.cpp
)
set(TEHPC_DEPEND_FILES ${TEHPC_SRC} ${TEHPC_HEADER} PARENT_SCOPE)
add_library(tehpc ${TEHPC_SRC})
# if linking to other libraries, e.g. MPI
#target_link_libraries(tehpc ${ARMADILLO_LIBRARIES})
#include "Model.hpp"
#include <cmath>
// constructor
Model::Model() :
nfi(0),
ekin(0), epot(0), temp(0),
rx(0,0),ry(0,0),rz(0,0),
vx(0,0),vy(0,0),vz(0,0),
fx(0,0),fy(0,0),fz(0,0)
{}
// destructor
Model::~Model() {
if (this->erg.is_open()) {
this->erg.close();
}
if (this->traj.is_open()) {
this->traj.close();
}
}
// reads a line for an integer value
void Model::read_line(int &el, std::ifstream &fin){
char line[BLEN];
fin.getline(line,BLEN,'#');
el = std::stoi(line); // Converts string to integer
fin.getline(line,BLEN);
}
//reads a line for a floating point value (double)
void Model::read_line(double &el, std::ifstream &fin){
char line[BLEN];
fin.getline(line,BLEN,'#');
el = std::stof(line); // Converts string to floating point
fin.getline(line,BLEN);
}
// reads input file and populates parameter struct
void Model::read_input(const std::string & fname){
std::ifstream fin(fname); // file input stream
read_line(this->natoms,fin);
read_line(this->mass,fin);
this->mass = this->mvsq2e * this->mass;
read_line(this->epsilon,fin);
read_line(this->sigma,fin);
read_line(this->rcut,fin);
read_line(this->box,fin);
read_line(this->nsteps,fin);
read_line(this->dt,fin);
read_line(this->nprint,fin);
this->rx.resize(this->natoms);
this->ry.resize(this->natoms);
this->rz.resize(this->natoms);
this->vx.resize(this->natoms);
this->vy.resize(this->natoms);
this->vz.resize(this->natoms);
this->fx.resize(this->natoms);
this->fy.resize(this->natoms);
this->fz.resize(this->natoms);
}
// reads restart file
void Model::read_restart(const std::string & fname) {
std::fstream rest;
rest.open(fname, std::ios::in);
if (rest.is_open()){
for (int i=0; i<this->natoms; ++i) {
rest >> this->rx[i];
rest >> this->ry[i];
rest >> this->rz[i];
}
for (int i=0; i<this->natoms; ++i) {
rest >> this->vx[i];
rest >> this->vy[i];
rest >> this->vz[i];
}
rest.close();
}
}
// init model
void Model::init() {
this->force();
this->compute_ekin();
}
// corrects distance between particles (x) by mapping it
// to the simulation box following periodic boundary conditions (pbc)
double Model::pbc(double x, const double boxby2){
while (x > boxby2) x -= 2.0*boxby2;
while (x < -boxby2) x += 2.0*boxby2;
return x;
}
// computes kinetic energy and temperature for ideal gas
// (equipartition theorem). Ekin = 3/2 Kboltz*T = 1/2 mv**2
void Model::compute_ekin() {
this->ekin=0.0;
for (int i=0; i<this->natoms; ++i) {
this->ekin += 0.5*this->mass*(this->vx[i]*this->vx[i] + this->vy[i]*this->vy[i] + this->vz[i]*this->vz[i]);
}
this->temp = 2.0*this->ekin/(3.0*this->natoms-3.0)/kboltz;
}
// computes force acting on every particle
void Model::force() {
double rcut2 = this->rcut * this->rcut; // cutoff radius squared
// zero energy and forces
this->epot=0.0;
std::fill(this->fx.begin(),this->fx.end(), 0);
std::fill(this->fy.begin(),this->fy.end(), 0);
std::fill(this->fz.begin(),this->fz.end(), 0);
for(int i=0; i < (this->natoms); ++i) {
for(int j=i; j < (this->natoms); ++j) {
// particles have no interactions with themselves
if (i==j) continue;
// get distance between particle i and j
double rx=pbc(this->rx[i] - this->rx[j], 0.5*this->box);
double ry=pbc(this->ry[i] - this->ry[j], 0.5*this->box);
double rz=pbc(this->rz[i] - this->rz[j], 0.5*this->box);
double r2 = rx*rx + ry*ry + rz*rz; // distance-squared
// compute force and energy only if particles are within cutoff
// distance
if (r2 < rcut2) {
double r = sqrt(r2); // distance between particles
// lennard-jones force
double ffac = -4.0*this->epsilon*(-12.0*pow(this->sigma/r,12.0)/r
+6*pow(this->sigma/r,6.0)/r);
// Lennard - Jones potential
this->epot += 4.0*this->epsilon*(pow(this->sigma/r,12.0)
-pow(this->sigma/r,6.0));
// split force into components
// Use Newton's 2nd Law (Fij = -Fji)
this->fx[i] += rx/r*ffac; this->fx[j] -= rx/r*ffac;
this->fy[i] += ry/r*ffac; this->fy[j] -= ry/r*ffac;
this->fz[i] += rz/r*ffac; this->fz[j] -= rz/r*ffac;
}
}
}
}
// integration algorithm: velocity verlet
// second-order central difference method
void Model::velverlet() {
// predictor stage: propagate velocities by half and positions by full step
for (int i=0; i<this->natoms; ++i) {
this->vx[i] += 0.5*this->dt * this->fx[i] / this->mass;
this->vy[i] += 0.5*this->dt * this->fy[i] / this->mass;
this->vz[i] += 0.5*this->dt * this->fz[i] / this->mass;
this->rx[i] += this->dt*this->vx[i];
this->ry[i] += this->dt*this->vy[i];
this->rz[i] += this->dt*this->vz[i];
}
// compute forces and potential energy
this->force();
// corrector stage: propagate velocities by another half step with new force
for (int i=0; i<this->natoms; ++i) {
this->vx[i] += 0.5*this->dt * this->fx[i] / this->mass;
this->vy[i] += 0.5*this->dt * this->fy[i] / this->mass;
this->vz[i] += 0.5*this->dt * this->fz[i] / this->mass;
}
}
// init dump
void Model::init_dump(const std::string & path,
const std::string & fname) {
std::string path_to_file = path + "/" + fname;
this->erg.open(path_to_file+".dat",std::fstream::out);
this->traj.open(path_to_file+".xyz",std::fstream::out);
}
// appends data to output.
void Model::dump() {
if ((this->nfi % this->nprint) == 0) {
// temperature and energy
if (this->erg.is_open()){
this->erg << this->nfi <<" "<<this->temp<<" "<<this->ekin<<" "<<this->epot<<" "<<this->ekin+this->epot<<"\n";
}
// atom positions
if (this->traj.is_open()){
this->traj << this->natoms <<" nfi="<<this->nfi<<" etot="<<this->ekin+this->epot<<"\n";
for (int i=0; i<this->natoms; ++i) {
this->traj << "Ar "<< this->rx[i]<< " "<< this->ry[i]<<" "<< this->rz[i]<< "\n";
}
}
}
}
#ifndef MODELHEADERDEF
#define MODELHEADERDEF
#include <fstream> // file streaming for input file
#include <vector>
// generic file- or pathname buffer length
#define BLEN 200
class Model
{
public:
// constructor and destructor
Model();
~Model();
// init
void init();
// input reading
void read_input(const std::string&);
void read_restart(const std::string&);
// model modification
double pbc(double , const double );
void force();
void velverlet();
void increment_time_step() { this->nfi++; }
// model output