1 #include "lmp2.h" 2 unwrap_molecules(struct NewAtomCoordinates * coord,struct Sys * sysinfo)3void 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