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 4634d308 authored by dkammer's avatar dkammer
Browse files

add week 10 group 2

parent b1586980
build
*~
\ No newline at end of file
# Assignment
This week we look into how we can use existing code to make our code faster/better/simpler/etc. Typically, we can use external libraries for this, as we have done with armadillo in the direct-stiffness-method example. For the Molecular Dynamics code in this assignment, we will look into BLAS. Please follow these steps:
1. read about BLAS:
1. BLAS: [https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms](https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms)
2. ATLAS: [https://en.wikipedia.org/wiki/Automatically_Tuned_Linear_Algebra_Software](https://en.wikipedia.org/wiki/Automatically_Tuned_Linear_Algebra_Software)
3. In case, the difference is note entirely clear: [https://stackoverflow.com/questions/17858104/what-is-the-relation-between-blas-lapack-and-atlas](https://stackoverflow.com/questions/17858104/what-is-the-relation-between-blas-lapack-and-atlas)
4. summarize on slides
2. However, using BLAS in your C++ code may not be straight-forward. For this reason, we will use the GSL library:
1. read about GSL:
1. [https://en.wikipedia.org/wiki/GNU_Scientific_Library](https://en.wikipedia.org/wiki/GNU_Scientific_Library)
2. [https://www.gnu.org/software/gsl/](https://www.gnu.org/software/gsl/)
2. summarize on slides
3. Link your code to BLAS and GSL:
1. in your `code/CMakeLists.txt` find BLAS and GSL:
1. `find_package(BLAS REQUIRED)`
2. `include_directories(${BLAS_INCLUDE_PATH})`
3. `find_package(GSL REQUIRED)`
4. `include_directories(${GSL_INCLUDE_DIRS})`
2. in your `code/src/CMakeLists.txt` use `${GSL_LIBRARIES} ${BLAS_LIBRARIES}` to link library to.
3. compile to check it works so far.
4. We want to use GSL in the Model to multiply a vector with a constant and add it to another vector (e.g. in `velverlet()` function).
1. Have a look at [https://www.gnu.org/software/gsl/doc/html/blas.html#level-1](https://www.gnu.org/software/gsl/doc/html/blas.html#level-1), if there is an appropriate function for this.
5. This is IMPORTANT: While we want to use GSL in our code, it is good practice to minimize the parts of our code that depend on an external library (so that we could change it easily). For this reason, we implemented a new class `PVector`, which is the only part of our code that depends on GSL (and includes GSL).
1. study the class `PVector`
2. summarize the purpose of each member and function on slides
3. change `code/src/CMakeLists.txt` to include the `PVector` in the compilation
6. Let's modify our code:
1. In `Model.hpp` change `rz`, `vz`, and `fz` to `PVector` instead of `std::vector<double>`. Important: do not change the x and y components, because we want to compare the performance.
2. Try to compile and fix all compilation errors
3. Run code and verify if still same result
4. In `velverlet()`, there is a block marked by "CHANGE THIS". Modify it to use the power of BLAS and the `PVector` class
5. Compile and run the code
6. Verify that the result is still the same
7. Note the time measured for computing the Y-component vs. the Z-component (printed to the terminal)
8. summarize on slides
7. For your discussion: what would be a good strategy to account for the possibility that someone does not have access to GSL? How would you change the code to provide an option to use or not use GSL/BLAS?
\ No newline at end of file
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})
# add BLAS and GSL here
# 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()
# ------------------------------------------------------
# test
option(TEHPC_TEST "tests" OFF)
if(${TEHPC_TEST})
enable_testing()
add_subdirectory(tests)
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) {
if (i%100 == 0)
std::cout << "step " << i << "/" << nb_increments << std::flush << std::endl;
my_model.increment_time_step();
my_model.velverlet();
my_model.compute_ekin();
my_model.dump();
}
std::cout << "end time stepping" << std::endl;
return 0;
}
#!/usr/bin/env python3
import sys
import numpy as np
import matplotlib.pyplot as plt
def plot(sname):
data = []
with open(sname+".dat",'r') as fl:
lines = fl.readlines()
for line in lines:
data.append([float(i) for i in line.strip().split()])
data = np.array(data)
fig = plt.figure()
ax = fig.add_subplot(311)
ax.plot(data[:,0],data[:,2],label='kin')
ax.legend(loc='best')
ax.set_ylabel('energy')
ax = fig.add_subplot(312)
ax.plot(data[:,0],data[:,3],label='pot')
ax.legend(loc='best')
ax.set_ylabel('energy')
ax = fig.add_subplot(313)
ax.plot(data[:,0],data[:,4],label='tot')
ax.legend(loc='best')
ax.set_xlabel('iter')
ax.set_ylabel('energy')
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(data[:,0],data[:,1],label='temperature')
ax.legend(loc='best')
ax.set_xlabel('iter')
ax.set_ylabel('temperature')
plt.show()
if __name__ == "__main__":
if len(sys.argv) < 2:
sys.exit('Missing argument! usage: ./basic_example.plot.py <sim-name>')
sname = str(sys.argv[1])
plot(sname)
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
# PVector.hpp
Model.hpp
)
set(TEHPC_SRC
# PVector.cpp
Model.cpp
)
set(TEHPC_DEPEND_FILES ${TEHPC_SRC} ${TEHPC_HEADER} PARENT_SCOPE)
add_library(tehpc ${TEHPC_SRC})
# if linking to other libraries, e.g. MPI
#include "Model.hpp"
#include <cmath>
#include <iostream>
#include <chrono>
// 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),
t1(0), t2(0) // to measure times
{}
// destructor
Model::~Model() {
// print measured times when done
std::cout << "t1 = " << this->t1 << std::endl;
std::cout << "t2 = " << this->t2 << std::endl;
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->vx.resize(this->natoms);
this->fx.resize(this->natoms);
this->ry.resize(this->natoms);
this->vy.resize(this->natoms);
this->fy.resize(this->natoms);
this->rz.resize(this->natoms);
this->vz.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;
}
void Model::force() {
this->epot = this->force(this->fx, this->fy, this->fz);
}
// computes force acting on every particle
double Model::force(std::vector<double> & fx,
std::vector<double> & fy,
std::vector<double> & fz) {
double rcut2 = this->rcut * this->rcut; // cutoff radius squared
// zero energy and forces
double epot = 0.0;
std::fill(fx.begin(),fx.end(), 0);
std::fill(fy.begin(),fy.end(), 0);
std::fill(fz.begin(),fz.end(), 0);