From: Ilya Chorny (ichorny_at_gmail.com)
Date: Mon Nov 12 2007 - 13:15:07 CST
Below is a Vega ZZ C script for calculating area/lipid. To use select
the lipid atoms for which you want to calulate the area per lipid
using the selection tool in Vega ZZ. You will need to know the number
of atoms per lipid. Edit the value of natoms accordingly to get the
proper normalization. Note the area per lipid in the presence of the
protein is a highly fluctuating quantity. I also do not image the box
so the values at the end of the box are not correct.
/*******************************************
**** VEGA ZZ C Script ****
**** Standard template ****
**** (c) 2005-2007, Alessandro Pedretti ****
*******************************************/
#define VGVLL
#include <vega.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
void Error(const char *Err)
{
VegaCmdEx("MessageBox \"%s.\" \"ERROR\" 16", Err);
}
int VllMain(VGP_VEGAINFO *VegaInfo)
{
ATOMO *Atm1;
ATOMO *Atm2;
int i, Ac, TrjTotFrm, j;
int xval, yval, binsize;
FILE *fp;
char buf[1024];
if (!*(VegaInfo -> BegAtm)) {
Error("No molecule");
return 0;
}
binsize = 6.5; //in angstroms
VegaCmd("Get TrjTotFrm");
TrjTotFrm = atoi(VegaInfo -> Result);
if (!TrjTotFrm) {
Error("Trajectory not open");
return 0;
}
if ((fp = fopen("pdens.txt", "w")) == NULL){
Error("Can't open file");
return 0;
}
//Determine MaxMin positions of the phosphates in the first frame
VegaCmdEx("TrjSel %d", 1);
int maxX, minX, maxY, minY;
maxX = 0;
minX = 0;
maxY = 0;
minY = 0;
for(Atm1 = *VegaInfo -> BegAtm; Atm1; Atm1 = Atm1 -> Ptr) {
if (Atm1 -> Active) {
xval = (int) Atm1 -> x;
if(xval > maxX)
maxX = xval;
if(xval < minX)
minX = xval;
yval = (int) Atm1 -> y;
if(yval > maxY)
maxY = yval;
if(yval < minY)
minY = yval;
}
}
// provide a buffer on either side of the min max values;
maxY += 10;
maxX += 10;
minX -= 10;
minY -= 10;
// allocate memory for distribution array
int lengthX, lengthY;
lengthX = (maxX - minX)/binsize+1;
lengthY = (maxY - minY)/binsize+1;
double **density;
density = calloc(lengthX, sizeof(double));
for(i=0;i<lengthX;i++)
density[i] = (double*) calloc(lengthY, sizeof(double));
// calculate distribution
int numframes, nleaflets, natoms, initframe;
initframe = 1;
numframes = TrjTotFrm;
nleaflets = 2;
natoms = 125;
for(j = initframe; j <= numframes; j+=1) {
VegaCmdEx("TrjSel %d", j);
for(Atm1 = *VegaInfo -> BegAtm; Atm1; Atm1 = Atm1 -> Ptr) {
if (Atm1 -> Active){
xval = ((int) Atm1 -> x - minX)/binsize;
yval = ((int) Atm1 -> y - minY)/binsize;
density[xval][yval] += 1.0;
// fprintf(fp, "%d %d\n", xval, yval);
}
}
}
for(i=0;i<lengthX;i++){
for(j=0;j<lengthY;j++){
if(density[i][j] > 0.0)
fprintf(fp, "%f\t",
((numframes-initframe-1.0)*binsize*binsize*nleaflets*natoms)/density[i][j]);
else
fprintf(fp, "%f\t", density[i][j]);
}
fprintf(fp, "\n");
}
fclose(fp);
return 0;
}
-- Ilya Chorny Ph.D.
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:45:31 CST