379 
379 
real arcdist;

380 
380 
gmx_rmpbc_t gpbc=NULL;

381 
381 


382 
// BEGIN CN MOD 


383 
// variables for modification by CN


384 
rvec ab,ac,bc,e,ce,eg,ef,g,f,cg,cf;


385 
float dab,def,deg,gdot,fdot,edot;


386 
// END CN MOD 


387 


388 

382 
389 
/* PBC added for centerofmass vector*/

383 
390 
/* Initiate the pbc structure */

384 
391 
memset(&pbc,0,sizeof(pbc));

...  ...  
501 
508 
direction[0],direction[1],direction[2]);*/

502 
509 
}

503 
510 

504 

if (bUnsat) {

505 

/* Using convention for unsaturated carbons */

506 

/* first get Sz, the vector from Cn to Cn+1 */

507 

rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);

508 

length = norm(dist);

509 

check_length(length, a[index[i]+j], a[index[i+1]+j]);

510 

svmul(1/length, dist, Sz);

511 


512 

/* this is actually the cosine of the angle between the double bond

513 

and axis, because Sz is normalized and the two other components of

514 

the axis on the bilayer are zero */

515 

if (use_unitvector)

516 

{

517 

sdbangle += acos(iprod(direction,Sz)); /*this can probably be optimized*/

518 

}

519 

else

520 

sdbangle += acos(Sz[axis]);

521 

} else {

522 

/* get vector dist(Cn1,Cn+1) for tail atoms */

523 

rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i1]+j]], dist);

524 

length = norm(dist); /* determine distance between two atoms */

525 

check_length(length, a[index[i1]+j], a[index[i+1]+j]);

526 


527 

svmul(1/length, dist, Sz);

528 

/* Sz is now the molecular axis Sz, normalized and all that */

529 

}

530 


531 

/* now get Sx. Sx is normal to the plane of Cn1, Cn and Cn+1 so

532 

we can use the outer product of Cn1>Cn and Cn+1>Cn, I hope */

533 

rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);

534 

rvec_sub(x1[a[index[i1]+j]], x1[a[index[i]+j]], tmp2);

535 

cprod(tmp1, tmp2, Sx);

536 

svmul(1/norm(Sx), Sx, Sx);

537 


538 

/* now we can get Sy from the outer product of Sx and Sz */

539 

cprod(Sz, Sx, Sy);

540 

svmul(1/norm(Sy), Sy, Sy);

541 


542 

/* the square of cosine of the angle between dist and the axis.

543 

Using the innerproduct, but two of the three elements are zero

544 

Determine the sum of the orderparameter of all atoms in group

545 

*/

546 

if (use_unitvector)

547 

{

548 

cossum[XX] = sqr(iprod(Sx,direction)); /* this is allowed, since Sa is normalized */

549 

cossum[YY] = sqr(iprod(Sy,direction));

550 

cossum[ZZ] = sqr(iprod(Sz,direction));

551 

}

552 

else

553 

{

554 

cossum[XX] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized */

555 

cossum[YY] = sqr(Sy[axis]);

556 

cossum[ZZ] = sqr(Sz[axis]);


511 
// BEGIN CN MOD 


512 
if (bUnsat) {


513 
//exactly as for the case without unsat, but here we select the ce vector to represent the hydrogen


514 
// place the 2 imaginary hydrogen atoms by enscribing a tetrahedron in a cube


515 
// a,b,c are the 3 backbone atoms, where a and b are in two cube corners and c is at the center


516 
rvec_sub(x1[a[index[i1]+j]], x1[a[index[i+1]+j]], ab);


517 
rvec_sub(x1[a[index[i]+j]], x1[a[index[i+1]+j]], ac);


518 
rvec_sub(x1[a[index[i]+j]], x1[a[index[i1]+j]], bc);


519 
// e is a point in the middle of the face opposite from the ab face


520 
for(m=0;m<DIM;m++){


521 
ce[m]=(ac[m]+bc[m])/2;

557 
522 
}


523 
rvec_add(x1[a[index[i]+j]],ce,e);


524 
//UNCOMMENT THE NEXT LINE TO PRINT OUT THE POSITION OF THE PSEUDOATOM FOR TESTING <<<<<<<<<< LOOK LOOK LOOK LOOK


525 
fprintf(stdout,"ATOM E %f %f %f\n",e[0],e[1],e[2]);


526 
//now take the angle between ce and z and this gets plugged into S_CD=1/2<3cos^2(ang)1>


527 
edot=ce[2]/norm(ce);


528 


529 
// using frameorder[ZZ] to store the result


530 
frameorder[ZZ] += 0.5 * (3 * edot*edot  1);


531 


532 
} else {


533 
// place the 2 imaginary hydrogen atoms by enscribing a tetrahedron in a cube


534 
// a,b,c are the 3 backbone atoms, where a and b are in two cube corners and c is at the center


535 
rvec_sub(x1[a[index[i1]+j]], x1[a[index[i+1]+j]], ab);


536 
rvec_sub(x1[a[index[i]+j]], x1[a[index[i+1]+j]], ac);


537 
rvec_sub(x1[a[index[i]+j]], x1[a[index[i1]+j]], bc);


538 
// e is a point in the middle of the face opposite from the ab face


539 
for(m=0;m<DIM;m++){


540 
ce[m]=(ac[m]+bc[m])/2;


541 
}


542 
rvec_add(x1[a[index[i]+j]],ce,e);


543 
dab=norm(ab);


544 
//the eg vector is the cross product of the ac and ce vectors (and ef is cross prod of bc and ce)


545 
cprod(ac,ce,eg);


546 
cprod(bc,ce,ef);


547 
deg=norm(eg);


548 
def=norm(ef);


549 
//progress along eg by the correct distance from e to find g (and along ef from e to find f)


550 
for(m=0;m<DIM;m++){


551 
g[m]=e[m]eg[m]*(dab/deg)/2;


552 
f[m]=e[m]ef[m]*(dab/def)/2;


553 
}


554 
//now find the cg and cf vectors, which are the CD vectors for S_CD


555 
rvec_sub(g,x1[a[index[i]+j]],cg);


556 
rvec_sub(f,x1[a[index[i]+j]],cf);


557 
//UNCOMMENT THE NEXT TWO LINES TO PRINT OUT THE POSITION OF THE PSEUDOATOM FOR TESTING <<<<<<<<<< LOOK LOOK LOOK LOOK


558 
fprintf(stdout,"ATOM G %f %f %f\n",g[0],g[1],g[2]);


559 
fprintf(stdout,"ATOM F %f %f %f\n",f[0],f[1],f[2]);


560 


561 
//now take the angle between cg and z  and also for cf and z  and this gets plugged into S_CD=1/2<3cos^2(ang)1>


562 
gdot=cg[2]/norm(cg);


563 
fdot=cf[2]/norm(cf);


564 


565 
// using frameorder[ZZ] to store the result  mult by 0.25 instead of 0.5 because we are adding two values


566 
frameorder[ZZ] += 0.25 * (3 * gdot*gdot  1);


567 
frameorder[ZZ] += 0.25 * (3 * fdot*fdot  1);


568 
}

558 
569 

559 

for (m = 0; m < DIM; m++)

560 

frameorder[m] += 0.5 * (3 * cossum[m]  1);

561 


562 

if (bSliced) {

563 

/* get average coordinate in box length for slicing,

564 

determine which slice atom is in, increase count for that

565 

slice. slFrameorder and slOrder are reals, not

566 

rvecs. Only the component [axis] of the order tensor is

567 

kept, until I find it necessary to know the others too

568 

*/

569 


570 

z1 = x1[a[index[i1]+j]][axis];

571 

z2 = x1[a[index[i+1]+j]][axis];

572 

z_ave = 0.5 * (z1 + z2);

573 

if (z_ave < 0)

574 

z_ave += box[axis][axis];

575 

if (z_ave > box[axis][axis])

576 

z_ave = box[axis][axis];

577 


578 

slice = (int)(0.5 + (z_ave / (*slWidth)))  1;

579 

slCount[slice]++; /* determine slice, increase count */

580 


581 

slFrameorder[slice] += 0.5 * (3 * cossum[axis]  1);

582 

}

583 

else if (permolecule)

584 

{

585 

/* store permolecule order parameter

586 

* To just track singleaxis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction)  1);

587 

* following is for Scd order: */

588 

(*slOrder)[j][i] += 1* (0.3333 * (3 * cossum[XX]  1) + 0.3333 * 0.5 * (3 * cossum[YY]  1));

589 

}

590 

if (distcalc)

591 

{

592 

/* bin order parameter by arc distance from reference group*/

593 

arcdist=acos(iprod(dvec,direction));

594 

(*distvals)[j][i]+=arcdist;

595 

}


570 
// END CN MOD 

596 
571 
} /* end loop j, over all atoms in group */

597 
572 

598 
573 
for (m = 0; m < DIM; m++)

...  ...  
702 
677 
for (atom = 1; atom < ngrps  1; atom++) {

703 
678 
fprintf(ord,"%12d %12g %12g %12g\n", atom, order[atom][XX],

704 
679 
order[atom][YY], order[atom][ZZ]);

705 

fprintf(slOrd,"%12d %12g\n", atom, 1 * (0.6667 * order[atom][XX] +

706 

0.333 * order[atom][YY]));


680 
// BEGIN CN MOD 


681 
// will place the order parameter value directly in the ZZ index.


682 
// not sure why I need to take the negative value....


683 
fprintf(slOrd,"%12d %12g\n", atom, 1 * order[atom][ZZ]);


684 
// END CN MOD 

707 
685 
}

708 
686 

709 
687 
ffclose(ord);
