Vega ZZ Script for calculating area/lipid.

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:47:14 CST