API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Nov 27th, 2014
171
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 10.88 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <cmath>
  3. #include <random>
  4. #include <cstdlib>
  5. #include <ctime>
  6. #include <new>
  7. using namespace std;
  8.  
  9. #include "mathlink.h"
  10. extern void mchshdpre6linux(float *wr0, int wr0len, float *wr1, int wr1len, float *wr2, int wr2len, float *plist, int plistlen, float solFactor, int poplen, int vectlen, int seed, int mcmax, int start, int stop);
  11. extern void randomsphere(float c[3], float r);
  12. extern int round_int(float r);
  13. extern int int_IntRange(int a, int b);
  14.  
  15.  
  16. void randomonsphere(float c[3], float r)
  17. {
  18. float pi = 3.14159265;
  19. float theta = ((float)rand() / (float)RAND_MAX) * 2 * pi;
  20. float u = 2 * ((float)rand() / (float)RAND_MAX) - 1;
  21. c[0] = r*cos(theta)*sqrt(1.0 - u*u);
  22. c[1] = r*sin(theta)*sqrt(1.0 - u*u);
  23. c[2] = r*u ;
  24.  
  25. }
  26.  
  27. int round_int(float r)
  28. {
  29. return (r > 0.0) ? (r + 0.5) : (r - 0.5);
  30. }
  31.  
  32. int rand_IntRange(int a, int b)
  33. {
  34. return rand() % (b - a + 1) + a;
  35. }
  36.  
  37. void mchshdpre6linux(float *wr0, int wr0len, float *wr1, int wr1len, float *wr2, int wr2len, float *plist, int plistlen, float solFactor, int poplen, int vectlen, int seed, int mcmax, int start, int stop)
  38. {
  39. /* This is the intial configuration for the two chains*/
  40.  
  41. float sigmaa = plist[0];
  42. float sigmab = plist[1];
  43. float bl = plist[2];
  44. float t = plist[3];
  45. //float bl = 3.93;
  46. //float t =42.6
  47. //float sigma = 3.93;
  48. float chain[4][3] = { { bl, 0.0, 1.0 }, { 0.0, 0.0, 0.0 }, { bl, 0.0, 0.0 }, { 0.0, 0.0, 0.0 } };
  49. float delr = 0.1;
  50. int ngrid = 2048;
  51. int temp[2][2] = { { 0, 0 }, { 0, 0 } };
  52. float sphere = 1.0;
  53. int rmin = start;
  54. int rmax = stop;
  55. float sr;
  56. float x1[3] = { 0., 0., 0. };
  57. float x2[3] = { 0., 0., 0. };
  58. float dir[3] = { 0., 0., 0. };
  59. float vr[3];
  60. float dis[3];
  61. float energy;
  62. int rIndex;
  63. float rFactor;
  64. float solE;
  65. float ljE;
  66. float gr[2048][3] = { 0.0, 0.0, 0.0 };
  67. float kb = 0.0019833794749;
  68. float entemp = 0.0;
  69. int dim[2];
  70. dim[0] = 2048;
  71. dim[1] = 3;
  72. //float hr[2048] = { 0 };
  73. float sigma[2][2]={ { 0, 0 }, { 0, 0 } };
  74. int choice1;
  75. int choice2;
  76. int choice3;
  77. float sited[3];
  78.  
  79. sigma[0][0] = 0.5*(sigmaa + sigmaa);
  80. sigma[0][1] = 0.5*(sigmaa + sigmab);
  81. sigma[1][0] = 0.5*(sigmab + sigmaa);
  82. sigma[1][1] = 0.5*(sigmab + sigmab);
  83.  
  84. srand(seed);
  85.  
  86. //const int numrows = poplen;
  87. int numcols = 3;
  88. //generate new array for vectors
  89. float **vectorlist = new float*[vectlen];
  90. for (int i = 0; i < vectlen; ++i)
  91. {
  92. vectorlist[i] = new float[3];
  93. }
  94. //generate new array orientations of dimers
  95. float **chainpop = new float*[poplen];
  96. for (int i = 0; i < poplen; ++i)
  97. {
  98. chainpop[i] = new float[3];
  99. }
  100.  
  101. for (int p = 0; p < poplen; p++)
  102. {
  103. randomonsphere(sited, bl);
  104. for (int c = 0; c < numcols; c++)
  105. {
  106. chainpop[p][c] = sited[c];
  107. }
  108. }
  109. for (int p = 0; p < vectlen; p++)
  110. {
  111. randomonsphere(sited, 1);
  112. for (int c = 0; c < numcols; c++)
  113. {
  114. vectorlist[p][c] = sited[c];
  115. }
  116. }
  117.  
  118. for (int i = 1; i <= mcmax; i++)
  119. {
  120. for (int n = rmin; n <= rmax; n++)
  121. {
  122. for (int j = 0; j <= 1; j++)
  123. {
  124. for (int k = 0; k <= 1; k++)
  125. {
  126. temp[k][j] = 0.0;
  127. }
  128. }
  129. /*randomonsphere(x1, bl);
  130. randomonsphere(x2, bl);*/
  131. choice1 = rand_IntRange(0, poplen-1);
  132. choice2 = rand_IntRange(0, poplen-1);
  133. choice3 = rand_IntRange(0, vectlen - 1);
  134. for (int q = 0; q < 3;q++)
  135. {
  136. x1[q] = chainpop[choice1][q];
  137. x2[q] = chainpop[choice2][q];
  138. dis[q] = vectorlist[choice3][q]*delr*n;
  139. }
  140. solE = 0.0;
  141. ljE = 0.0;
  142. int da = rand() % 2;
  143. int db = rand() % 2;
  144. if (da == 1)
  145. {
  146. for (int j = 0; j <= 2; j++)
  147. {
  148. chain[0][j] = 0.0;
  149. chain[1][j] = x1[j];
  150. }
  151. }
  152. else if (da == 0)
  153. {
  154. for (int j = 0; j <= 2; j++)
  155. {
  156. chain[1][j] = 0.0;
  157. chain[0][j] = x1[j];
  158. }
  159. }
  160. if (db == 1)
  161. {
  162. for (int j = 0; j <= 2; j++)
  163. {
  164. chain[2][j] = 0.0;
  165. chain[3][j] = x2[j];
  166. }
  167. }
  168. else if (db ==0)
  169. {
  170. for (int j = 0; j <= 2; j++)
  171. {
  172. chain[3][j] = 0.0;
  173. chain[2][j] = x2[j];
  174. }
  175. }
  176. /*dis[0] = dir[0]*delr*n;
  177. dis[1] = dir[1]*delr*n;
  178. dis[2] = dir[2]*delr*n;*/
  179. energy = 0.0;
  180.  
  181. for (int j = 0; j <= 1; j++)
  182. {
  183. for (int k = 0; k <= 1; k++)
  184. {
  185. for (int s = 0; s <= 2; s++)
  186. {
  187. vr[s] = chain[k + 2][s] - chain[j][s] + dis[s];
  188. sr = 0.0;
  189. }
  190. for (int s = 0; s <= 2; s++)
  191. {
  192. sr = vr[s] * vr[s] + sr;
  193. }
  194. sr = sqrt(sr);
  195. temp[j][k] = round_int(sr / delr) - 1;
  196. rIndex = round_int(sr / delr) - 1;
  197. rFactor = sr / delr - rIndex - 1;
  198. if (rIndex == 0)
  199. {
  200. if (j == 0 && k == 0){
  201. solE = solE + solFactor*wr0[0];
  202. }
  203. else if (j == 0 && k == 1){
  204. solE = solE + solFactor*wr1[0];
  205. }
  206. else if (j == 1 && k == 0){
  207. solE = solE + solFactor*wr1[0];
  208. }
  209. else if (j == 1 && k == 1){
  210. solE = solE + solFactor*wr2[0];
  211. }
  212. }
  213. else if (rIndex < wr0len)
  214. {
  215. if (j == 0 && k == 0)
  216. {
  217. solE = solE + solFactor*(wr0[rIndex] * (1.0 - rFactor) + rFactor*wr0[rIndex + 1]);
  218. }
  219. else if (j == 0 && k == 1)
  220. {
  221. solE = solE + solFactor*(wr1[rIndex] * (1.0 - rFactor) + rFactor*wr1[rIndex + 1]);
  222. }
  223. else if (j == 1 && k == 0)
  224. {
  225. solE = solE + solFactor*(wr1[rIndex] * (1.0 - rFactor) + rFactor*wr1[rIndex + 1]);
  226. }
  227. else if (j == 1 && k == 1)
  228. {
  229. solE = solE + solFactor*(wr2[rIndex] * (1.0 - rFactor) + rFactor*wr2[rIndex + 1]);
  230. }
  231.  
  232. }
  233. //if (sr*sr >= sigma*sigma)
  234. if (sr*sr >= sigma[j][k]*sigma[j][k])
  235. {
  236. entemp = 0.0;
  237. }
  238. else
  239. {
  240. entemp = 1000000000.0;
  241. }
  242. energy = energy + entemp;
  243.  
  244. }
  245. }
  246.  
  247. for (int j = 0; j <= 1; j++)
  248. {
  249. for (int k = 0; k <= 1; k++)
  250. {
  251. if (j == 0 && k == 0){
  252. gr[temp[j][k]][0] = gr[temp[j][k]][0] + exp(-(energy + solE) / (kb*t))*n*n;
  253. }
  254. else if (j == 0 && k == 1){
  255. gr[temp[j][k]][1] = gr[temp[j][k]][1] + exp(-(energy + solE) / (kb*t))*n*n;
  256. }
  257. else if (j == 1 && k == 0){
  258. gr[temp[j][k]][1] = gr[temp[j][k]][1] + exp(-(energy + solE) / (kb*t))*n*n;
  259. }
  260. else if (j == 1 && k == 1){
  261. gr[temp[j][k]][2] = gr[temp[j][k]][2] + exp(-(energy + solE) / (kb*t))*n*n;
  262. }
  263.  
  264. }
  265. }
  266.  
  267. }
  268. }
  269.  
  270.  
  271. // free dynamic memory
  272. for(int i = 0; i < vectlen; i++)
  273. {
  274. delete[] vectorlist[i];
  275. }
  276. for(int i = 0; i < poplen; i++)
  277. {
  278. delete[] chainpop[i];
  279. }
  280. delete[] vectorlist;
  281. delete[] chainpop;
  282. MLPutReal32Array(stdlink, (float *)gr, (int *)dim,(char **)0, 2);
  283.  
  284. return;
  285. }
  286.  
  287. #if WINDOWS_MATHLINK
  288.  
  289. #if __BORLANDC__
  290. #pragma argsused
  291. #endif
  292.  
  293. int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
  294. {
  295. char buff[512];
  296. char FAR * buff_start = buff;
  297. char FAR * argv[32];
  298. char FAR * FAR * argv_end = argv + 32;
  299.  
  300. hinstPrevious = hinstPrevious; /* suppress warning */
  301.  
  302. if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;
  303. MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
  304. return MLMain( (int)(argv_end - argv), argv);
  305. }
  306.  
  307. #else
  308.  
  309. int main(int argc, char* argv[])
  310. {
  311. return MLMain(argc, argv);
  312. }
  313.  
  314. #endif
  315.  
  316. :Begin:
  317. :Function: mchshdpre6linux
  318. :Pattern: mchshdpre6linux[a1_List, a2_List, a3_List, p1_List, b_Real, c_Integer, d_Integer, e_Integer, f_Integer, g_Integer, h_Integer]
  319. :Arguments: {a1, a2, a3, p1, b, c, d, e, f, g, h}
  320. :ArgumentTypes: {Real32List, Real32List, Real32List, Real32List, Real32, Integer, Integer, Integer, Integer, Integer, Integer}
  321. :ReturnType: Manual
  322. :End:
  323.  
  324. # This makefile can be used to build all or some of the sample
  325. # programs. To build all of them, use the command
  326. # 'make all'. To build one, say addtwo, use the command
  327. # 'make addtwo'.
  328.  
  329. # Portions of this makefile require the use of GNU make.
  330. # see http://www.gnu.org/software/make for more information.
  331.  
  332. VERSION=10.0
  333. MLINKDIR = $(shell pwd)/..
  334. SYS = Linux-x86-64
  335. CADDSDIR = $//usr/local/Wolfram/Mathematica/10.0/SystemFiles/Links/MathLink/DeveloperKit/Linux-x86-64/CompilerAdditions
  336. EXTRA_CFLAGS=-m64
  337. CXXFLAGS =-std=c++0x
  338.  
  339. INCDIR = ${CADDSDIR}
  340. LIBDIR = ${CADDSDIR}
  341.  
  342. MPREP = ${CADDSDIR}/mprep
  343. RM = rm
  344.  
  345. CC = /usr/bin/cc
  346. CXX = /usr/bin/c++
  347.  
  348. BINARIES = mchshdpre6linux
  349.  
  350. all : $(BINARIES)
  351.  
  352. mchshdpre6linux : twositetwochaindimertm.o twositetwochaindimer.o
  353. ${CXX} ${CXXFLAGS} ${EXTRA_CFLAGS} -I${INCDIR} twositetwochaindimertm.o twositetwochaindimer.o -L${LIBDIR} -lML64i4 -lm -lpthread -lrt -lstdc++ -ldl -luuid -o $@
  354.  
  355.  
  356. .c.o :
  357. ${CC} -c ${EXTRA_CFLAGS} -I${INCDIR} $<
  358.  
  359. twositetwochaindimertm.c : twositetwochaindimer.tm
  360. ${MPREP} $? -o $@
  361.  
  362.  
  363.  
  364. clean :
  365. @ ${RM} -rf *.o *tm.c $(BINARIES)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!