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

week 11

parent 09f80d69
build
*~
\ No newline at end of file
# Assignment
This week we look a bit more into data structures. We work on the same code we have used over the past couple of weeks. The objective here is to learn about iterators and use them to loop over built-in data structures. For this assignment, follow these steps:
1. read about iterators:
1. [https://www.geeksforgeeks.org/introduction-iterators-c/](https://www.geeksforgeeks.org/introduction-iterators-c/)
2. [https://www.geeksforgeeks.org/difference-between-iterators-and-pointers-in-c-c-with-examples/](https://www.geeksforgeeks.org/difference-between-iterators-and-pointers-in-c-c-with-examples/)
3. [https://www.cplusplus.com/reference/iterator/](https://www.cplusplus.com/reference/iterator/)
2. summarize what you learnt on 1-2 slides
3. evaluate how you could use iterators in the `velverlet` function and summarize on 1-2 slides
4. make the changes and summarize problems/insights/results on slides (keep a copie of for the moment)
5. measure the time for the `velverlet` and compare to previous code. Helpful code to measure time:
1. `#include <chrono> // to time function execution`
2. `auto t1 = std::chrono::high_resolution_clock::now(); // get current time`
3. `auto dt = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count() // time passed between t1 and t2`
As a reminder: The code is provided in the folder `code`. It can be compile via `cmake`. Follow these steps:
1. `cd code`
2. `mkdir build`
3. `cd build`
4. `ccmake ..`
5. configure and generate
6. `make`
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 )
# OPENMP
#option(TEHPC_USE_OPENMP "Use OPENMP: velverlet on 3 threads: 1 dimension per thread (shared memory)" OFF)
#if (TEHPC_USE_OPENMP)
# set(CMAKE_CXX_FLAGS "-std=c++11 -fopenmp" CACHE STRING "" FORCE)
#endif (TEHPC_USE_OPENMP)
# config file: provide cmake variables to cpp code
#configure_file(tehpc_config.hpp.in
# "${CMAKE_CURRENT_BINARY_DIR}/src/tehpc_config.hpp" @ONLY)
# include directories
set(TEHPC_INCLUDE_DIRS
"${CMAKE_CURRENT_BINARY_DIR}/src"
"${CMAKE_CURRENT_SOURCE_DIR}/src"
)
include_directories(${TEHPC_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()
# ------------------------------------------------------
# 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
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 "tehpc_config.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), t3(0), t4(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;
std::cout << "t3 = " << this->t3 << std::endl;
std::cout << "t4 = " << this->t4 << 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->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;
}
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);
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;