1 #include "lmp2.h"
2 
unwrap_molecules(struct NewAtomCoordinates * coord,struct Sys * sysinfo)3 void unwrap_molecules(struct NewAtomCoordinates *coord,struct Sys *sysinfo)
4 {
5   int iflag[5],nflag[5],icase;
6   register int i,iax,imol;
7   int min_true,max_true;
8 
9   if (trueflag) {
10 
11     /* Use trueflags to "unwrap" molecules */
12 
13     for (imol=0; imol < sysinfo->no_molecules; imol++) {
14       for (iax=0; iax < 3; iax++) {
15 
16 	min_true =  1000;
17 	max_true = -1000;
18 	for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
19 	     i++) {
20 	  if (coord[i].truef[iax] < min_true) min_true = coord[i].truef[iax];
21 	  if (coord[i].truef[iax] > max_true) max_true = coord[i].truef[iax];
22 	}
23 
24 	if ((min_true > 0) || (max_true < 0)) {
25 	  for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
26 	       i++) coord[i].truef[iax] -= min_true;
27 	}
28 
29       } /* end loop over iax */
30     } /* end loop over imol */
31   }
32   else {
33 
34     /* Use coordinates to "unwrap" molecules */
35 
36     for (imol=0; imol < sysinfo->no_molecules; imol++) {
37 
38       for (iax=0; iax < 3; iax++) {
39 
40 	for (i=0; i<5; i++) { iflag[i] = 0; nflag[i] = 0;}
41 
42 	for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
43 	     i++) {
44 	  if (coord[i].fract[iax] <= 0.20)
45 	    nflag[0]++;
46 	  else if ((coord[i].fract[iax] > 0.20) && (coord[i].fract[iax] <= 0.40))
47 	    nflag[1]++;
48 	  else if ((coord[i].fract[iax] > 0.40) && (coord[i].fract[iax] <= 0.60))
49 	    nflag[2]++;
50 	  else if ((coord[i].fract[iax] > 0.60) && (coord[i].fract[iax] <= 0.80))
51 	    nflag[3]++;
52 	  else if ((coord[i].fract[iax] > 0.80) && (coord[i].fract[iax] <= 1.00))
53 	    nflag[4]++;
54 	}
55 
56 	for (i=0; i<5; i++) { if (nflag[i] > 0) iflag[i] = 1; }
57 
58 	icase = 10000*iflag[0] + 1000*iflag[1] + 100*iflag[2] + 10*iflag[3] + iflag[4];
59 
60 	switch (icase) {
61 	 case 10001:
62 	   if (nflag[0] > nflag[4]) {
63 	     for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
64 		  i++) if (coord[i].fract[iax] > 0.80) coord[i].fract[iax] -= 1.0;
65 	   }
66 	   else {
67 	     for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
68 		  i++) if (coord[i].fract[iax] <= 0.20) coord[i].fract[iax] += 1.0;
69 	   }
70 	   break;
71          case 11001:
72 	   for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
73 		i++) if (coord[i].fract[iax] > 0.80) coord[i].fract[iax] -= 1.0;
74 	   break;
75          case 10011:
76 	   for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
77 		i++) if (coord[i].fract[iax] < 0.20) coord[i].fract[iax] += 1.0;
78 	   break;
79          case 11101:
80 	   for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
81 		i++) if (coord[i].fract[iax] > 0.80) coord[i].fract[iax] -= 1.0;
82 	   break;
83          case 10111:
84 	   for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
85 		i++) if (coord[i].fract[iax] < 0.20) coord[i].fract[iax] += 1.0;
86 	   break;
87          case 11011:
88 	   if ((nflag[0]+nflag[1]) > (nflag[3]+nflag[4])) {
89 	     for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
90 		  i++) if (coord[i].fract[iax] > 0.60) coord[i].fract[iax] -= 1.0;
91 	   }
92 	   else {
93 	     for (i=sysinfo->molinfo[imol].start; i < sysinfo->molinfo[imol].end;
94 		  i++) if (coord[i].fract[iax] <= 0.40) coord[i].fract[iax] += 1.0;
95 	   }
96 	   break;
97          default:
98 	   break;
99 	} /* end switch */
100 
101       } /* end loop over imol*/
102     } /* end loop over iax */
103 
104   } /* end if on trueflag */
105 }
106