1 /* -*- tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
2 vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
3
4 MDAnalysis --- https://www.mdanalysis.org
5 Copyright (c) 2006-2016 The MDAnalysis Development Team and contributors
6 (see the file AUTHORS for the full list of names)
7
8 Released under the GNU Public Licence, v2 or any higher version
9
10 Please cite your use of MDAnalysis in published work:
11
12 R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
13 D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
14 MDAnalysis: A Python package for the rapid analysis of molecular dynamics
15 simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
16 Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
17
18 N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
19 MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
20 J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26 #include <math.h>
27 #include <time.h>
28 #include <sys/types.h>
29
30 #ifdef __unix__
31 #include <unistd.h>
32 #endif
33
34 #ifdef _WIN32
35 #include <io.h>
36 #endif
37
38 #define EPSILON 1e-8
39
40 typedef struct {
41 int index;
42 double value;
43 } IVWrapper;
44
trmIndex(int row,int col)45 inline int trmIndex(int row, int col) { // array index for triangular matrix
46 return (row) > (col) ? ((row)+1)*(row)/2+(col) : ((col)+1)*(col)/2+(row);
47 }
48
printarray(double * arr,int len)49 void printarray(double* arr, int len) {
50 for (int i=0; i<len; i++) {
51 printf("%.4f ",arr[i]); }
52 printf("\n");
53 }
54
cmp_ivwrapper(const void * ptra,const void * ptrb)55 int cmp_ivwrapper(const void* ptra, const void* ptrb) {
56 IVWrapper* a = (IVWrapper*) ptra;
57 IVWrapper* b = (IVWrapper*) ptrb;
58 if (a->value == b->value)
59 return(0);
60 else {
61 if (a->value < b->value)
62 return(-1);
63 else
64 return(1);
65 }
66 }
67
68 // Euclidean distance
ed(double * d_coords,int i,int j,int dim)69 double ed(double* d_coords, int i, int j, int dim) {
70 double d = 0.0;
71 double delta = 0.0;
72 int idxi = i*dim;
73 int idxj = j*dim;
74 int idxik = 0;
75 int idxjk = 0;
76
77 for (int k=0; k<dim; k++) {
78 idxik = idxi+k;
79 idxjk = idxj+k;
80 delta = d_coords[idxik] - d_coords[idxjk];
81 d = d + delta*delta;
82 }
83 return(sqrt(d));
84 }
85
stress(double * s,double * d_coords,int dim,int elemsn)86 double stress(double *s, double *d_coords, int dim, int elemsn) {
87 double denom = 0.0;
88 double numer = 0.0;
89 double dab = 0.0;
90 double delta = 0.0;
91 int k = 0;
92
93 for (int i=0; i<elemsn; i++) {
94 for (int j=0; j<i; j++) {
95 dab = ed(d_coords, i, j, dim);
96 denom += s[k];
97 delta = dab - s[k];
98 numer += delta*delta / s[k];
99 //printf("%d, dab %.3f, denom %.3f, numer %.3f, s[k]: %.3f\n", k, dab, denom, numer,s[k]);
100 k++;
101 //printf("a %d, b %d, dab %.3f, rab %.3f, delta %.3f\n",i,j,dab,s[k],dab-s[k]);
102 }
103 k++;
104 }
105 return( numer/denom );
106 }
107
neighbours(double * s,int nelem,double cutoff,int * s_indeces,int * ioffsets,int * js)108 int neighbours(double *s, int nelem, double cutoff, int* s_indeces, int* ioffsets, int* js) {
109 int idx = 0;
110 int global_counter = 0;
111 int offset_counter = 0;
112 ioffsets[0] = 0;
113 for (int i=0;i<nelem;i++) {
114 offset_counter = 0;
115 for (int j=0;j<nelem;j++) {
116 idx = trmIndex(i,j);
117 if (s[idx] < cutoff) {
118 s_indeces[global_counter] = idx; //XXX: replace numbers with addresses
119 js[global_counter] = j;
120 offset_counter++;
121 global_counter++;
122 }
123 }
124 ioffsets[i+1] = ioffsets[i] + offset_counter;
125 }
126 return global_counter;
127 }
128
nearest_neighbours(double * s,int nelem,int kn)129 int* nearest_neighbours(double *s, int nelem, int kn) {
130 IVWrapper* ivpairs = (IVWrapper*) malloc((nelem-1)*sizeof(IVWrapper));
131 int* neighbours = (int*) malloc(nelem*kn*sizeof(int));
132 int totk = 0;
133 int k = 0;
134 for (int i=0;i<nelem;i++) {
135 k = 0;
136 for (int j=0;j<nelem;j++) {
137 if (i!=j) {
138 ivpairs[k].index = j;
139 ivpairs[k].value = s[trmIndex(i,j)];
140 k++;
141 }
142 }
143 qsort(ivpairs, nelem-1, sizeof(IVWrapper), cmp_ivwrapper);
144 for (int k=0;k<kn;k++) {
145 neighbours[totk] = ivpairs[k].index;
146 totk++;
147 }
148 }
149 free(ivpairs);
150 return neighbours;
151 }
152
neighbours_stress(double * s,double * d_coords,int dim,int elemsn,double rco)153 double neighbours_stress(double *s, double *d_coords, int dim, int elemsn, double rco) {
154 double denom = 0.0;
155 double numer = 0.0;
156 double dab = 0.0;
157 double delta = 0.0;
158 int k = 0;
159
160 for (int i=0; i<elemsn; i++) {
161 for (int j=0; j<i; j++) {
162 dab = ed(d_coords, i, j, dim);
163 if (s[k] <= rco || dab < s[k]) {
164 denom += s[k];
165 delta = dab - s[k];
166 numer += delta*delta / s[k];
167 }
168 //printf("%d, dab %.3f, denom %.3f, numer %.3f, s[k]: %.3f\n", k, dab, denom, numer,s[k]);
169 k++;
170 //printf("a %d, b %d, dab %.3f, rab %.3f, delta %.3f\n",i,j,dab,s[k],dab-s[k]);
171 }
172 k++;
173 }
174 return( numer/denom );
175 }
176
177
CStochasticProximityEmbedding(double * s,double * d_coords,double rco,int nelem,int dim,double maxlam,double minlam,int ncycle,int nstep,int stressfreq)178 double CStochasticProximityEmbedding(
179 double* s,
180 double* d_coords,
181 double rco,
182 int nelem,
183 int dim,
184 double maxlam,
185 double minlam,
186 int ncycle,
187 int nstep,
188 int stressfreq) {
189
190 int a = 0, b = 0, idxa = 0, idxb = 0, idxak = 0, idxbk = 0;
191 double dab = 0.0, rab = 0.0;
192 double lam = maxlam;
193 double t = 0.0;
194 double finalstress = 0.0;
195
196 /* random init of d */
197 srand(time(NULL)+getpid()*getpid());
198 for (int i=0; i<nelem*dim; i++) {
199 d_coords[i] = (double) rand() / (double) RAND_MAX;
200 }
201
202 /* start self organization */
203 for (int i=0; i<ncycle; i++) {
204 for (int j=0; j<nstep; j++) {
205
206 a = rand() % nelem;
207 while(1) {
208 b = rand() % nelem;
209 if (b != a) break;
210 }
211
212 dab = ed(d_coords, a, b, dim);
213 rab = s[trmIndex(a, b)];
214
215 if (rab <= rco || (rab > rco && dab < rab)) {
216 idxa = a * dim;
217 idxb = b * dim;
218 t = lam * 0.5 * (rab - dab) / (dab + EPSILON);
219
220 for (int k=0; k<dim; k++) {
221 idxak = idxa+k;
222 idxbk = idxb+k;
223 d_coords[idxak] = d_coords[idxak] + t*(d_coords[idxak] - d_coords[idxbk]);
224 d_coords[idxbk] = d_coords[idxbk] + t*(d_coords[idxbk] - d_coords[idxak]);
225 }
226 }
227 }
228 lam = lam - (maxlam - minlam) / (double)(ncycle - 1);
229 //if (i % stressfreq == 0 && i != 0 && stressfreq > 0)
230 //printf("Cycle %d - Residual stress: %.3f, lambda %.3f\n", i, neighbours_stress(s, d_coords, dim, nelem, rco),lam);
231 }
232 finalstress = neighbours_stress(s, d_coords, dim, nelem, rco);
233 //printf("Calculation finished. - Residual stress: %.3f\n", finalstress);
234 return(finalstress);
235 }
236