Project

General

Profile

g_order.patch

The patch listed in the mailing list post - Chris Neale, 02/27/2013 10:35 PM

View differences:

gmx_order.c 2011-06-09 17:57:11.586591000 -0400
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 center-of-mass 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(Cn-1,Cn+1) for tail atoms */
523
	  rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
524
	  length = norm(dist);      /* determine distance between two atoms */
525
	  check_length(length, a[index[i-1]+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 Cn-1, Cn and Cn+1 so
532
	   we can use the outer product of Cn-1->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[i-1]+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[i-1]+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[i-1]+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 PSEUDO-ATOM 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[i-1]+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[i-1]+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 PSEUDO-ATOM 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[i-1]+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 per-molecule order parameter
586
		 *  To just track single-axis 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);