main.c

#include "functions.c"
#define nequil 20000

int main(int argc, char **argv) {
  gsl_error_handler_t * old_handler=gsl_set_error_handler_off();

  //variable init
  int s,n;
  double dx, dy, dz, fLx, fRx, Ro, kbt;
  double fLprev=0.0, fRprev=0.0, fxi=0.0;
  double sumsq, sumhelix, rms, xL, zL, xR, zR;
  int accepted = 0;
  int accepted_chain = 0;
  int mcstep = 0;
  double De;
  double totalE;
  char filename[64];
  int i; int ranN;

  Int n_rxz = 50, n_y = 20;
  const Doub rxzs[] = {0.00011248105,0.0002445754,0.00039970507,0.0005818891,0.00079584776,0.0010471249,0.0013422316,0.0016888149,0.0020958562,0.0025739038,0.0031353464,0.0037947343,0.0045691562,0.0054786826,0.0065468857,0.0078014503,0.0092748916,0.011005398,0.01303782,0.015424834,0.018228308,0.021520902,0.025387954,0.029929694,0.035263835,0.041528631,0.048886457,0.057528017,0.067677289,0.079597329,0.093597089,0.11003942,0.1293505,0.15203084,0.17866828,0.20995325,0.24669659,0.28985067,0.34053399,0.40006019,0.46997212,0.55208183,0.64851749,0.76177862,0.8948008,1.0510318,1.234521,1.4500242,1.703127,2.0003895} ;
  const Doub ys[] = {0.0,0.026315789,0.052631579,0.078947368,0.10526316,0.13157895,0.15789474,0.18421053,0.21052632,0.23684211,0.26315789,0.28947368,0.31578947,0.34210526,0.36842105,0.39473684,0.42105263,0.44736842,0.47368421,0.5};
  const Doub zs[] = {4444.023768762637,18.88892548072377,9.38939720668971,6.224797428434006,4.64639979362448,3.704417707101945,3.081499104845088,2.641061619901676,2.314947155543283,2.065916002981771,1.872125268121191,1.719656318399189,1.598849512924905,1.502766524163137,1.426655550158602,1.367558164648733,1.323681400760863,1.293701660950426,1.276394087296795,1.270757915410431,2044.243278803969,18.88408158859048,9.38730399650497,6.224903938461017,4.647503320256135,3.705189116558555,3.081360028065836,2.640376679788338,2.31447243414313,2.066058036552988,1.872645066171977,1.719994954578368,1.598695976702739,1.502323123611393,1.426390235458616,1.367725643331337,1.324087303201154,1.293920533577516,1.276205226301456,1.270362660962849,1250.806411212508,18.88271006966062,9.38713221896558,6.224853119016275,4.647482023170942,3.705178214514184,3.081353589208329,2.640372494870193,2.314469600283751,2.066056099898398,1.872643697824008,1.719993889769956,1.598695057003728,1.502322314659383,1.426389570074546,1.367725131710047,1.324086887762504,1.293920122878848,1.27620476701084,1.270362173305297,859.154338651823,18.88025808458066,9.3868255012118,6.224762154181137,4.647443579396997,3.70515847180665,3.081342109734204,2.640365214734729,2.31446467388215,2.066052591333905,1.872641091209715,1.719991881358994,1.598693457969278,1.502321001999788,1.426388460168162,1.367724165287276,1.324086020908443,1.293919321388276,1.276204002607205,1.270361420926604,628.1449329623229,18.87621838491461,9.38631998264162,6.224612219291265,4.647380211964219,3.705125929234885,3.081323187618503,2.640353214532619,2.314456553435305,2.06604680796578,1.872636794568698,1.719988570769321,1.598690822185326,1.502318838266682,1.426386630643989,1.367722572269228,1.324084592013545,1.29391800023674,1.276202742594693,1.270360180737231,477.3820299787004,18.86987722500922,9.38552597668789,6.224376693291318,4.647280666928689,3.705074806579771,3.081293461678468,2.640334362546388,2.314443796387019,2.0660377223931,1.87263004460625,1.719983369871507,1.598686681395225,1.502315439060359,1.426383756475714,1.367720069648347,1.324082347228845,1.293915924714947,1.276200763122302,1.270358232406055,372.3980281234154,18.86023430931338,9.38431739879626,6.224018129828177,4.647129110797076,3.704996970725647,3.081248202294875,2.64030565908759,2.314424372775837,2.066023888834248,1.872619767182027,1.719975451018978,1.598680376645026,1.502310263435416,1.426379380268279,1.367716259155673,1.32407892931582,1.293912764520255,1.276197749172167,1.270355265871205,295.9496659100354,18.8458950974918,9.38251766316237,6.223484039502605,4.6469033428159,3.704881016348003,3.081180776476953,2.640262897093081,2.314395435496547,2.066003279479286,1.872604455717019,1.719963653350057,1.59867098367447,1.502302552647524,1.426372860464289,1.367710582163021,1.324073837196381,1.293908056354792,1.276193258885367,1.270350846224396,238.450048169062,18.82492606401584,9.37988029443954,6.222701068115512,4.646572323921687,3.704710994374537,3.081081907717956,2.64020019242667,2.31435300232384,2.06597305796093,1.872582002862693,1.719946353067064,1.598657209624981,1.502291245352224,1.426363299635563,1.367702257241991,1.324066369941675,1.29390115213648,1.276186674168216,1.270344365094924,194.1415166565066,18.79466285327312,9.37606237116619,6.221566976398806,4.646092765726763,3.704464655146723,3.080938652852756,2.640109334392972,2.314291516135786,2.065929266012072,1.872549467636977,1.719921283971879,1.59863725013097,1.502274860274001,1.426349445264765,1.367690193763802,1.324055549269213,1.293891147336412,1.276177132344753,1.270334973376371,159.3560730174715,18.75146102930403,9.37058832584069,6.219939617259655,4.645404429261428,3.704111023883374,3.080732989063556,2.639978888391377,2.314203237139002,2.065866390340622,1.872502753498481,1.719885289337107,1.598608591717003,1.502251333960606,1.426329552538633,1.367672872415245,1.3240400123589,1.293876781871174,1.276163431631472,1.270321488187085,131.6455912288793,18.69038167073241,9.36280092193684,6.217621832861386,4.644423657657415,3.703607059206022,3.080439864559294,2.639792957220944,2.314077403567971,2.065776764347259,1.87243616357843,1.719833979045171,1.598567738723163,1.50221779655054,1.426301194707831,1.367648180045149,1.324017863695702,1.293856303104003,1.276143900471668,1.270302264258287,109.3134510716391,18.60481242071784,9.3517954645866,6.214340845240217,4.643034500765159,3.702893055579302,3.080024512723613,2.639529472911391,2.313899073645739,2.065649742322712,1.872341786604671,1.719761256115613,1.598509836144724,1.502170262033288,1.426261001106713,1.367613181494361,1.323986470368955,1.293827276554929,1.276116216991339,1.27027501622535,91.146865656007,18.48604369564307,9.33633263142966,6.209720287789452,4.641076573719024,3.701886330892635,3.079438757180788,2.639157843998466,2.313647529341509,2.065470560554632,1.872208649683738,1.719658663286053,1.598428149016002,1.50210320062619,1.426204295417381,1.367563804512183,1.323942179376884,1.293786324472798,1.276077159642322,1.270236573185171,76.25623300120948,18.32285655510958,9.31472507582338,6.203242538623448,4.638328518628466,3.700472583749423,3.078615938818427,2.638635719960396,2.313294078778612,2.065218768215086,1.872021550419812,1.719514481921912,1.598313344241979,1.502008948745862,1.426124596459792,1.367494404865522,1.32387992739661,1.293728764936315,1.276022262951548,1.27018253984354,63.97467929354123,18.10124065473142,9.28469384085065,6.194198533810951,4.634485609292662,3.698494104735689,3.077463969247702,2.637904550348551,2.312799035435776,2.064866068223458,1.871759448681018,1.719312490477512,1.59815250067209,1.50187689544193,1.426012929673597,1.367397166411185,1.323792702384298,1.293648113869054,1.275945342478663,1.270106828912613,53.79300271804739,17.80444969058214,9.24319387911956,6.1816219129538,4.629129672104652,3.695733779347882,3.075855852765951,2.636883505709387,2.312107574237465,2.064373351053782,1.87139325492486,1.719030255782773,1.597927745870939,1.501692361518141,1.425856878181295,1.367261274026104,1.323670801173712,1.29353539826575,1.275837839625044,1.270001016162168,45.31624730801022,17.41370069082058,9.18621567574355,6.164204504742807,4.62168915096885,3.691893535434523,3.07361681768892,2.635461185798762,2.311144062537083,2.063686628630928,1.870882793983114,1.718636784761194,1.597614381088226,1.501435056850942,1.425639276099305,1.367071774288291,1.323500806492418,1.293378209666771,1.275687918527289,1.269853451395152,38.23393637169188,16.90987502927567,9.10858514848257,6.140191409947884,4.6113869571384,3.686565597056035,3.070506962956311,2.633484365423847,2.309804333667659,2.06273147675914,1.870172646697588,1.718089303393773,1.597178306402399,1.501076960035537,1.425336411295256,1.366808007592001,1.323264178602612,1.293159400691005,1.275479222172758,1.269648033890623,32.29918489019944,16.2764713153507,9.00381013834273,6.107257351263711,4.597173996090365,3.679194690291959,3.06619808029894,2.630742824141786,2.307945210616333,2.06140546945779,1.869186472889653,1.717328850056178,1.596572492433961,1.500579409284551,1.424915557200168,1.366441453325975,1.322935318633145,1.292855291901001,1.275189160830649,1.269362527298218,27.3137327861345,15.50366225414573,8.86406071673353,6.06237320080546,4.577647473142658,3.66902937647184,3.060243116101077,2.626949093605388,2.305370403036471,2.059567936768028,1.867819295151774,1.716274268283651,1.595732159735602,1.499889120900599,1.424331590171365,1.365932773618893,1.322478909535063,1.292433208481415,1.27478656001785,1.268966243720359,23.11701763525701,14.59260411176132,8.68041706686451,6.001684867427139,4.550956131515425,3.655061294679436,3.052036720449283,2.621711788729216,2.301811713628641,2.057026206539541,1.865927075084187,1.714814057903703,1.594568219024821,1.498932758169551,1.423522366618033,1.365227767301536,1.321846272850451,1.291848102604424,1.274228434002738,1.268416866284537,19.57806050672284,13.55843470252218,8.44355165867717,5.920448007148243,4.514700808856704,3.635952859074716,3.040765746635172,2.614501180560518,2.296904353086612,2.053517324941946,1.863312752392961,1.71279539677062,1.592958394631579,1.497609557075565,1.422402429563322,1.364251850567604,1.320970391307561,1.29103793399526,1.273455570306107,1.267656099890306,16.58934783146831,12.43028772467472,8.14498620902358,5.813090776804804,4.465849777838219,3.609957800394146,3.025349862707732,2.604606041667513,2.290155188748711,2.048684148051844,1.859707772718433,1.710009493005403,1.590735306560177,1.495781377812687,1.420854486830756,1.362902562087583,1.319759135501826,1.289917374456353,1.272386506892909,1.266603737043495,14.06215389014037,11.24761564329901,7.778916531629116,5.67350532187948,4.400705508479791,3.574846121283415,3.004375343181616,2.591081942432059,2.280903155553285,2.04204477727894,1.854748069628985,1.706172307960461,1.587670644837693,1.493259404128828,1.418717958073762,1.361039449852898,1.318086097180599,1.288369270235286,1.27090935434927,1.265149597137656,11.92291944175057,10.05380217732017,7.344300519983593,5.49567551830159,4.314981524819037,3.527856787006847,2.976030636293323,2.572693729815437,2.268272213911859,2.032954797658169,1.847943624167268,1.700899684934676,1.583454483629995,1.489786587008405,1.415773744262253,1.358470546044814,1.315778278401531,1.286233142738077,1.268870762375207,1.263142644793563,10.11041581884286,8.88931923914522,6.846548086857071,5.274696812598814,4.204067331817951,3.46571784459356,2.938061509907518,2.547859751824355,2.251119737088746,2.020563050237846,1.838641345320198,1.69367618521931,1.57766885453869,1.485014892774945,1.411724276211944,1.354934505406424,1.312599736434314,1.283289848811722,1.266061165900239,1.260376427443514,8.57350091599081,7.786634371003915,6.297991818466018,5.008098087815526,4.063557441279364,3.384790991317388,2.887776721488502,2.51461339524599,2.227987837965437,2.003764304926223,1.825982498137262,1.683817746962116,1.569755118304352,1.478476600357949,1.406167941557975,1.350077439391934,1.308230189223794,1.279241396536022,1.262195305111666,1.256569826457315,7.269326569441334,6.767908509512578,5.716622138673369,4.697147826064439,3.890068444557068,3.28140482680516,2.822149753598181,2.47060883623607,2.197072637296263,1.98115692129559,1.808858831697612,1.670430077715076,1.558975760832316,1.469549532416717,1.398567346073223,1.343423654611207,1.302237648941393,1.273684906391369,1.256886970765564,1.251342061595615,6.161893891682929,5.845213566663496,5.123309001817539,4.347631674111181,3.682243542411447,3.152417681012512,2.738069989828225,2.413210285277469,2.156236736816169,1.951020944790714,1.785876579407025,1.652368049649149,1.544373466761286,1.457417574541729,1.388211829799451,1.334340134203264,1.294044563398099,1.266079984364455,1.249617128861375,1.244181067785454,5.220879363726339,5.02229344461708,4.538449623275055,3.96961765960623,3.441670092648332,2.995974534521038,2.632781749964471,2.33971015084157,2.103100628687162,1.911341596090893,1.755344924511943,1.628206833823431,1.524734259778533,1.441030649229167,1.374176672709339,1.321996060443654,1.282888030201282,1.255709619360064,1.239695291632644,1.234405033481748,4.420673304361129,4.296865646921773,3.979164845455442,3.576083150377159,3.173314639210595,2.812293405114693,2.504492637678417,2.247707506899303,2.035252134933346,1.859911582399395,1.715314288207399,1.596243030536564,1.498567510775582,1.419072982136196,1.355285029444665,1.305321484293313,1.267776889165733,1.241636582922187,1.226215654466523,1.221118473544686,3.739586016414578,3.662788518951279,3.457740953332154,3.180801271855295,2.885153030255091,2.604190938713006,2.353029825852428,2.135627791637769,1.950596239608041,1.7945452304273,1.663695795930394,1.554551783099672,1.464123113424181,1.389955091585054,1.330084213067302,1.282973932576692,1.247452481545006,1.222660798799582,1.208012717280758,1.203167312718001,3.159187038492235,3.111782904903642,2.981367466887151,2.796214730032048,2.586983077415792,2.377044607374384,2.180313807879772,2.00327568014762,1.847817508605321,1.713415737766839,1.598487872026155,1.501128534747135,1.4194722273109,1.35185208855055,1.296855067894598,1.253328409999177,1.220365999100144,1.197288836908278,1.183626507414714,1.179102755435029,2.663750659645279,2.634648660472319,2.552777932065938,2.431937244024078,2.288788744445496,2.138069199764931,1.990391481961125,1.852221463626865,1.726846207698245,1.615474076717038,1.518107694809417,1.434132410426157,1.362676160735429,1.302814414852911,1.253679384673626,1.2145136383158,1.18469313623724,1.163734542292988,1.151295339765777,1.147171512259482,2.23978709393338,2.22204262050393,2.171294892125281,2.094136032066082,1.999210669980175,1.895092147008268,1.788899450056235,1.685803102926959,1.589145615394356,1.50083371901095,1.421769468798009,1.352212832043223,1.292048871089762,1.240969416051519,1.198589928540772,1.164521496255334,1.138413513621155,1.11997801925075,1.109002973620738,1.105359089905754,1.875644254453187,1.864925254909756,1.833899276369054,1.785665725109689,1.72455170297047,1.655242951657057,1.582078493578706,1.508643348543409,1.437641075864384,1.370952971531573,1.309786507353022,1.254844279493207,1.206476740273385,1.164804949804829,1.129812352742319,1.10141009744049,1.079481840612851,1.063913519631177,1.054612455623893,1.05151894898415,1.561170633546147,1.554780843285187,1.536124528007738,1.506642369422629,1.46844270470602,1.4239666559039,1.375671928101061,1.325798306718577,1.276235131924771,1.228478179535019,1.183648376081099,1.142544186456022,1.105706150525747,1.073480153759303,1.046072666679796,1.023595598765825,1.006100872580335,0.993605930290506,0.986111680789079,0.98361426406349,1.287435934861231,1.283698790735416,1.2727197441301,1.255163220413118,1.232036295162569,1.204565328501669,1.174066209330197,1.141831462530756,1.109047740669015,1.076746986556972,1.045787146204599,1.016854712661946,0.99048088807984,0.967064454940501,0.946896387776814,0.93018308692046,0.917066554150034,0.907640798073366,0.901964323448905,0.90006883009697,1.046512805685817,1.044385750824601,1.038109739629603,1.027989419461613,1.014498647101212,0.998236508219277,0.979877294032122,0.960121910929447,0.939656415751567,0.919120843199053,0.899089063377163,0.880058609923712,0.862448418211731,0.846602105819443,0.832794589147861,0.821240230888294,0.812101178491648,0.805494978774751,0.801500891182008,0.800164563455265,0.831328962489023,0.830163608334859,0.826714979622518,0.821121813796184,0.813603981317327,0.804447146919905,0.793984291686445,0.7825762897400327,0.7705935011729823,0.7583998366123223,0.7463401296026925,0.7347310650278482,0.7238554494887177,0.7139593104029183,0.7052511668673592,0.697902793326972,0.6920508549681959,0.6877988935793777,0.6852192561973548,0.684354668855369,0.6356006639625196,0.6349946572914095,0.633197812024478,0.6302724736813106,0.6263187471068436,0.6214692045397268,0.6158822328943179,0.6097346239522851,0.6032140074344772,0.5965116487687285,0.5898160056394096,0.58330728873925,0.5771531281027031,0.5715053251890847,0.5664975820469672,0.5622440443335125,0.5588384712328333,0.55635384614913,0.5548422604837018,0.5543349326051079,0.4538502851386648,0.4535562997619108,0.4526835554258818,0.451259293654035,0.449327633811652,0.4469477484290555,0.4441914957513148,0.4411406685852029,0.4378840296491638,0.4345142975197579,0.4311252269680016,0.4278088975561414,0.4246532898823213,0.4217401944911034,0.4191434678903499,0.4169276257363867,0.4151467462179602,0.4138436470847846,0.4130492969800624,0.4127824246553239,0.2814884477645783,0.2813580928918876,0.2809708387281658,0.2803379837189325,0.2794779226598059,0.2784155170053175,0.2771812594724776,0.2758102740965782,0.274341198262527,0.2728149949976703,0.2712737422269739,0.2697594413457658,0.2683128811600887,0.2669725858649992,0.2657738681265249,0.2647480012260338,0.2639215181551171,0.2633156408467645,0.2629458395256539,0.2628215204113488,0.1149094727941127,0.114857859819115,0.1147044720116266,0.1144536183542596,0.1141123343387267,0.1136901687719663,0.1131988969603254,0.112652171047912,0.1120651202076885,0.1114539145989415,0.1108353075056737,0.1102261699081673,0.1096430310118016,0.1091016370915314,0.1086165395497591,0.1082007214700726,0.1078652703040792,0.107619102752032,0.107468746456937,0.1074181818593529,-0.04847242075283665,-0.04849020145954791,-0.0485430530176775,-0.04862951772877249,-0.04874721227219698,-0.04889289542062745,-0.04906255981209834,-0.04925154490655481,-0.04945466764693438,-0.04966636687590149,-0.04988085724031325,-0.05009228814836378,-0.0502949033232357,-0.0504831966095459,-0.05065205991776888,-0.05079691951733261,-0.05091385729067963,-0.05099971401875702,-0.05105217226598677,-0.05106981695581261,-0.2103036709601284,-0.2103088499919367,-0.2103242453079567,-0.2103494354926346,-0.210383731149224,-0.2104261939517087,-0.2104756625471122,-0.210530784568218,-0.2105900538471743,-0.2106518517796664,-0.2107144916810371,-0.2107762649023163,-0.2108354874366864,-0.2108905457453279,-0.2109399405646176,-0.210982327522123,-0.2110165534838485,-0.2110416876763027,-0.2110570467703714,-0.2110622132757164,-0.3714896977391157,-0.3714909331167248,-0.3714946055203799,-0.3715006146860003,-0.3715087965593517,-0.3715189277861822,-0.3715307318236167,-0.3715438865040588,-0.3715580328430511,-0.371572784848671,-0.3715877400630232,-0.3716024905469128,-0.3716166340073477,-0.3716297847643418,-0.3716415842586037,-0.3716517108148605,-0.3716598883963664,-0.3716658941139573,-0.3716695642870463,-0.3716707988932697,-0.5324489082287029,-0.5324491409886932,-0.5324498329183776,-0.5324509651402463,-0.5324525067648712,-0.5324544157340715,-0.532456639968874,-0.5324591187908806,-0.5324617845781827,-0.5324645646105721,-0.5324673830536444,-0.5324701630276203,-0.5324728287044253,-0.5324753073758252,-0.5324775314362322,-0.5324794402261464,-0.5324809816860233,-0.5324821137755348,-0.5324828056195962,-0.5324830383499758,-0.6933394835362483,-0.6933395167632748,-0.6933396155379813,-0.6933397771659747,-0.6933399972383466,-0.6933402697519489,-0.693340587273161,-0.6933409411406762,-0.6933413217017764,-0.6933417185756481,-0.6933421209365545,-0.6933425178091391,-0.6933428983678044,-0.6933432522320009,-0.6933435697493701,-0.693343842259022,-0.6933440623277636,-0.6933442239528405,-0.6933443227256602,-0.6933443559520343} ;

  VecDoub lek_rxz(n_rxz, rxzs);
  VecDoub lek_y(n_y, ys);
  MatDoub lek_z(n_rxz,n_y, zs);
  Bilin_interp lekner(lek_rxz, lek_y, lek_z);
  
  //default parameters
  int N             = 16;
  int Nl            = 32;
  double ci_charge  = 1.0;
  double ep         = 1.0;
  double h          = 1.0;
  double htheta     = 1.0;
  double Lmax       = 0.8;
  int T             = 200000;
  double kf         = 1.0;
  double R          = 2.0;
  int posOut        = 1000;
  int dataOut      = 1000;
  int seed          = 1;
  double step       = 0.03;
  double stepChain  = 0.003;
  double turns      = 0.0;
  double amp        = 0.0;
  double wv         = 0.0;
  double Ly         = 24.0;
  int kc            = 0;
  double sigma      = 0.30;
  double sigma_c    = 0.15;
  double ebp        = 0.0;

  //import simulation parameters
  for(int i = 1; i < argc;i++) {
    if ( strcmp(argv[i], "-N") == 0 )
      N  = atoi(argv[++i]);
    else if ( strcmp(argv[i], "-Nl") == 0)
      Nl = atoi(argv[++i]);
    else if ( strcmp(argv[i], "-q") == 0 )
      ci_charge  = atof(argv[++i]);
    else if ( strcmp(argv[i], "-e") == 0 )
      ep = atof(argv[++i]);
    else if ( strcmp(argv[i], "-h") == 0 )
      h  = atof(argv[++i]);
    else if ( strcmp(argv[i], "-hth") == 0 )
      htheta = atof(argv[++i]);
    else if ( strcmp(argv[i], "-L") == 0 )
      Lmax   = atof(argv[++i]);
    else if ( strcmp(argv[i], "-T") == 0 )
      T      = atoi(argv[++i]);
    else if ( strcmp(argv[i], "-kf") == 0 )
      kf     = atof(argv[++i]);
    else if ( strcmp(argv[i], "-R") == 0 )
      R      = atof(argv[++i]);
    else if ( strcmp(argv[i], "-po") == 0 )
      posOut = atoi(argv[++i]);
    else if ( strcmp(argv[i], "-fo") == 0 )
      dataOut = atoi(argv[++i]);
    else if ( strcmp(argv[i], "-s") == 0 )
      seed     = atoi(argv[++i]);
    else if ( strcmp(argv[i], "-D") == 0 )
      step     = atof(argv[++i]);
    else if ( strcmp(argv[i], "-Dc") == 0 )
      stepChain= atof(argv[++i]);
    else if ( strcmp(argv[i], "-tr") == 0 )
      turns    = atof(argv[++i]);
    else if ( strcmp(argv[i], "-f") == 0 )
      strcpy(filename,argv[++i]); 
    else if ( strcmp(argv[i], "-A") == 0 )
      amp      = atof(argv[++i]);
    else if ( strcmp(argv[i], "-wv") == 0 )
      wv       = atof(argv[++i]);
    else if ( strcmp(argv[i], "-Ly") == 0 )
      Ly       = atof(argv[++i]);
    else if ( strcmp(argv[i], "-kc") == 0 )
      kc       = atoi(argv[++i]);
    else if ( strcmp(argv[i], "-sig") == 0 )
      sigma    = atof(argv[++i]);
    else if ( strcmp(argv[i], "-sigc") == 0 )
      sigma_c  = atof(argv[++i]);
    else if ( strcmp(argv[i], "-ebp") == 0 )
      ebp     = atof(argv[++i]);
  }
 
  //inits
  double qL = -0.5*N*ci_charge/(1.0*Nl);
  if(step == 0.0) {
    step       = 0.00075+0.0084*exp((kf-0.02842)/0.0326);
    stepChain  = step*0.01;
  }

  //init arrays
  double *q; 
  double **x; double **xn;
  q = (double *)malloc((N+2*Nl)*sizeof(double));
  x = (double **)malloc((N+2*Nl)*sizeof(double *));
  xn = (double **)malloc((N+2*Nl)*sizeof(double *));

  for(int i=0;i<N+2*Nl;i++)
    {
      x[i] = (double *)malloc(3*sizeof(double));
      xn[i] = (double *)malloc(3*sizeof(double));
      if(x[i] == NULL || xn[i] == NULL)
	{ printf("Out of memory\n"); }
    }

  //set random seed
  settable(1234567890987654321ULL, seed*1000 + 123456123456123456ULL, 362436362436362436ULL, 1066149217761810ULL);


  //open required output files
  char fname [64];
  FILE *pos;
  if (posOut != 0) {
    sprintf(fname, "%s.xyz",filename); pos=fopen(fname,"a");
  }
  char fnamedata [64];
  FILE *data;
  if (dataOut != 0) {
    sprintf(fnamedata, "%s.dat",filename); data=fopen(fnamedata,"a");
    //print headers
    fprintf(data,"%%seed\tt\tN\tNl\tci_charge\tep\th\ththeta\tLmax\tkf\tR\tturns\tfx\tRo\tRl\tA\twv\tA0\tstep\tstepC\tenergy\tsigma\tsigma_c\thelix\n");
    fflush(data);
  }



  //initialize perturbed chains
  if(kc != 0) {
    double phi1;
    double phi2;

    //linear portion
    for(i=N; i < N+2*Nl; i++) {
      if(i < N+Nl) {
        q[i]    = qL*0.1;
        x[i][0] = xn[i][0] = -0.5*R;
        x[i][1] = xn[i][1] = (i-N)*Ly/Nl;
        x[i][2] = xn[i][2] = 0.0;
      }
      else {
        q[i]    = qL*0.1;
        x[i][0] = xn[i][0] = 0.5*R;
        x[i][1] = xn[i][1] =  (i-N-Nl)*Ly/Nl;
        x[i][2] = xn[i][2] = 0.0;
      }
    }

    phi1 = twoPI*ran_u();
    phi2 = twoPI*ran_u();

    for(n=1; n <= kc; n++) {
      for(i=N; i < N+2*Nl; i++) {
        if(i < N+Nl) {
          x[i][0]  += amp*sin(twoPI*kc*(i-N)/Nl+phi1)/kc;
          xn[i][0] += amp*sin(twoPI*kc*(i-N)/Nl+phi1)/kc;
        }
        else {
          x[i][0]  += amp*sin(twoPI*kc*(i-N-Nl)/Nl+phi2)/kc;
          xn[i][0] += amp*sin(twoPI*kc*(i-N-Nl)/Nl+phi2)/kc;
        }
      }
    }
  }
  else if(turns != 0.0) { // Helix init
    for(i=N; i < N+2*Nl; i++) {
      if(i < N+Nl) {
        q[i]    = qL*0.1;
        x[i][0] = xn[i][0] = 0.5*R*sin(twoPI*turns*(i-N)/Nl);
        x[i][1] = xn[i][1] = (i-N)*Ly/Nl;
        x[i][2] = xn[i][2] = 0.5*R*cos(twoPI*turns*(i-N)/Nl);
      }
      else {
        q[i]    = qL*0.1;
        x[i][0] = xn[i][0] = 0.5*R*sin(twoPI*turns*(i-N-Nl)/Nl + PI);
        x[i][1] = xn[i][1] =  (i-N-Nl)*Ly/Nl;
        x[i][2] = xn[i][2] = 0.5*R*cos(twoPI*turns*(i-N-Nl)/Nl + PI);
      }
    }
  }
  else { // Wave perturbed init
    for(i=N; i < N+2*Nl; i++) {
      if(i < N+Nl) {
        q[i]    = qL*0.1;
        x[i][0] = xn[i][0] = -0.5*R+amp*cos(twoPI*(i-N)*wv/Nl);
        x[i][1] = xn[i][1] = (i-N)*Ly/Nl;
        x[i][2] = xn[i][2] = 0.0;
      }
      else {
        q[i]    = qL*0.1;
        x[i][0] = xn[i][0] = 0.5*R;
        x[i][1] = xn[i][1] =  (i-N-Nl)*Ly/Nl;
        x[i][2] = xn[i][2] = 0.0;
      }
    }
  }

  //initialize particles
  int w = 0;
  do {
    for(i=0; i<N;i++) {
      q[i]    = ci_charge*0.1;

      x[i][0] = xn[i][0] = ran_xz(0.5*R);
      x[i][1] = xn[i][1] = ran_y(Ly);
      x[i][2] = xn[i][2] = ran_xz(0.25*R);
    }
    w++;
  } while (w < 5000 && isfinite(energy(x,q,ep,sigma,sigma_c,h,htheta,Lmax,N,Nl,Ly,lekner)) == 0);

  //MC loop
  for(s = 0; s < T; s++) {

    if(posOut != 0 && s >= nequil && s % posOut == 0) {
      fprintf(pos,"%d\n",N+2*Nl);
      fprintf(pos,"N%d Nl%d q%f ep%f h%f Lmax%f T%d kbt%f R%f pos%d\n",N,Nl,ci_charge,ep,h,Lmax,T,kbt,R,posOut);
      for(int i=0;i<N+2*Nl;i++) {
        
        if(i >= N && i <N+Nl){
          fxi = 0.0;
          for(int j=0;j<N+2*Nl;j++){
            if(i == j || (j >= N && j <N+Nl)){
              continue;
            }
            fxi += lekner_fx(q[i],x[i][0],x[i][1],x[i][2],q[j],x[j][0],x[j][1],x[j][2],Ly);
          }
        } else {
          fxi = NAN;
        }
    
        fprintf(pos,"%d %10.10f %f %f %10.10f %f %d %d %d\n",on_chain(i,N)?10:1,x[i][0],x[i][1],x[i][2],fxi,R,(int)wv,seed,s);
      }

      fflush(pos);
    }

    if(dataOut != 0 && s >= nequil && s % dataOut == 0) {
      //initialize output vars
      fLx = 0.0;
      fRx = 0.0;
      xL  = 0.0;
      zL  = 0.0;
      xR  = 0.0;
      zR  = 0.0;
      sumsq= 0.0;
      sumhelix=0.0;

      //calculate total energy
      totalE = energy(x,q,ep,sigma,sigma_c,h,htheta,Lmax,N,Nl,Ly,lekner);

      //calculate forces
      for(int i=0; i<N+2*Nl; i++) {
        for(int j=0; j< N+2*Nl;j++) {
          if(i==j) { continue; }
          
          if(i >= N && i < N+Nl) {//if i is on left chain, add to fLx
            fLx += lekner_fx(q[i],x[i][0],x[i][1],x[i][2],q[j],x[j][0],x[j][1],x[j][2],Ly);
                //+  rep_fx(ep, x[i][0],x[i][1],x[i][2],x[j][0],x[j][1],x[j][2],Ly);
            fLx += linked(i,j,N,Nl) ? spring_fx(h, Lmax, x[i][0],x[i][1],x[i][2],x[j][0],x[j][1],x[j][2],Ly) : 0.0;
          }
          else if(i >= N+Nl) {//if i is on the right chain, add to fRx
            fRx += lekner_fx(q[i],x[i][0],x[i][1],x[i][2],q[j],x[j][0],x[j][1],x[j][2],Ly);
                //+  rep_fx(ep, x[i][0],x[i][1],x[i][2],x[j][0],x[j][1],x[j][2],Ly);
            fRx += linked(i,j,N,Nl) ? spring_fx(h, Lmax, x[i][0],x[i][1],x[i][2],x[j][0],x[j][1],x[j][2],Ly) : 0.0;
          }
        }
      }
 
      fLprev = fLx;
      fRprev = fRx;
      
      //calculate average position of lines
      for(int i=N; i<N+Nl; i++) {
        xL += x[i][0]/Nl;
        zL += x[i][2]/Nl;
      }
      for(int i=N+Nl; i<N+2*Nl; i++) {
        xR += x[i][0]/Nl;
        zR += x[i][2]/Nl;
      }
      Ro = sqrt((xR-xL)*(xR-xL)+(zR-zL)*(zR-zL));

      //calculate rms amplitude of left line
      for(int i=N; i<N+Nl; i++) {
        sumsq += (x[i][0]-xL)*(x[i][0]-xL);
      }
      rms = sqrt(2.0*sumsq/Nl);

      double currPhi, nextPhi;      
      //calculate total delta angle (helicity)
      for(int i=N; i<N+Nl-1; i++) {
        currPhi = atan2(x[i][2]-x[i+Nl][2], x[i][0]-x[i+Nl][0]);
        nextPhi = atan2(x[i+1][2]-x[i+1+Nl][2], x[i+1][0]-x[i+1+Nl][0]);
        sumhelix += nextPhi - currPhi;
      }

      fprintf(data,"%d\t%d\t%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", seed, s, N, Nl, ci_charge, ep, h, htheta, Lmax, kf, R, turns, 0.5*(fRx-fLx), Ro, xL, rms, wv,amp,step,stepChain,totalE,sigma,sigma_c,sumhelix);
      fflush(data);
    } //end dataOut
    if(s < nequil || stepChain == 0.0) {
      for(n = 0; n < N; n++) {
        kbt = kf;

        ranN = ran_particle(N);
        dx = ran_du(); dy = ran_du(); dz = ran_du();

        xn[ranN][0] += step*dx;
        xn[ranN][1] += step*dy;
        xn[ranN][2] += step*dz;

        De = delta_u(x,xn,q,ranN,ep,sigma,sigma_c,ebp,h,htheta,Lmax,N,Nl,Ly,lekner);

        if(De < 0 || exp(-De/kbt) > ran_u()) {
          if(s >= nequil) accepted += 1;
          x[ranN][0] = xn[ranN][0];
          x[ranN][1] = xn[ranN][1];
          x[ranN][2] = xn[ranN][2];
        }
        else {
          xn[ranN][0] = x[ranN][0];
          xn[ranN][1] = x[ranN][1];
          xn[ranN][2] = x[ranN][2];
        }

        if(s >= nequil) mcstep += 1;
      } //end particle only sweep
    }
    else {
      for(n = 0; n < N+2*Nl; n++) {
        kbt = kf;

        ranN = ran_particle(N+2*Nl);

        dx = ran_du(); dy = ran_du(); dz = ran_du();

        if(ranN < N) {
          xn[ranN][0] += step*dx;
          xn[ranN][1] += step*dy;
          xn[ranN][2] += step*dz;
        } else {
          xn[ranN][0] += stepChain*dx;
          xn[ranN][1] += is_endpoint(ranN,N,Nl) ? 0.0 : stepChain*dy;
          xn[ranN][2] += stepChain*dz;
        }

        De = delta_u(x,xn,q,ranN,ep,sigma,sigma_c,ebp,h,htheta,Lmax,N,Nl,Ly,lekner);

        if(De < 0 || exp(-De/kbt) > ran_u()) {
          if(s >= nequil) accepted += 1;
          if(ranN >= N) accepted_chain +=1;
          x[ranN][0] = xn[ranN][0];
          x[ranN][1] = xn[ranN][1];
          x[ranN][2] = xn[ranN][2];
        }
        else {
          xn[ranN][0] = x[ranN][0];
          xn[ranN][1] = x[ranN][1];
          xn[ranN][2] = x[ranN][2];
        }
        if(s >= nequil) mcstep += 1;
      } //end all sweep
    } //end particle/all if split

    if(s < nequil && s % 100 == 0) {
      for(i=0; i<N+2*Nl;i++) {
        if(i<N) {
          q[i] = ci_charge*(0.1+0.9*s/(1.0*(nequil-100)));
        }
        else {
          q[i] = qL*(0.1+0.9*s/(1.0*(nequil-100)));
        }
      }
    }
    
  } //end MC loops

  //close files
  if(posOut != 0) { fclose(pos); }
  if(dataOut != 0) { fclose(data); }

  for(int i=0; i< N+2*Nl; i++)
    {
      free(x[i]);
      free(xn[i]);
    }
  free(x); free(xn); free(q);

  printf("%f\t%f\n", (accepted-accepted_chain)/(mcstep*1.0*N/(N+2*Nl)), accepted_chain/(mcstep*2.0*Nl/(N+2*Nl)));

}