* And Hey:
* Gromacs Runs On Most of All Computer Systems
*/
-#ifndef _complex_h
-#define _complex_h
+#ifndef _gmxcomplex_h
+#define _gmxcomplex_h
#include <math.h>
#include "typedefs.h"
gmx_genpr.c gmx_eneconv.c gmx_vanhove.c gmx_wheel.c \
addconf.c addconf.h gmx_tune_pme.c gmx_membed.c \
calcpot.c calcpot.h edittop.c \
- linearsearch.c linearsearch.h dens_filter.c dens_filter.h \
+ dens_filter.c dens_filter.h \
powerspect.c powerspect.h interf.h \
gmx_densorder.c gmx_hydorder.c \
binsearch.h binsearch.c
-#include <types/simple.h>
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ * $Id: densorder.c,v 0.9
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.0
+ *
+ * Copyright (c) 1991-2001
+ * BIOSON Research Institute, Dept. of Biophysical Chemistry
+ * University of Groningen, The Netherlands
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
+ *
+ * And Hey:
+ * Gyas ROwers Mature At Cryogenic Speed
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <stdio.h>
+#include "types/simple.h"
/*Make range-array (Permutation identity) for sorting */
void rangeArray(int *ar,int size)
{
-int i;
+ int i;
for (i=0;i<size;i++){
ar[i]=i;
}
void Swap (real *v1, real *v2)
{
-real temp;
-temp = *v1;
-*v1 = *v2;
-*v2 = temp;
+ real temp;
+ temp = *v1;
+ *v1 = *v2;
+ *v2 = temp;
}
void insertionSort(real *arr, int *perm, int startndx, int endndx, int direction) {
- int i, j;
+ int i, j;
if(direction>=0){
- for (i = startndx; i <=endndx; i++) {
- j = i;
+ for (i = startndx; i <=endndx; i++) {
+ j = i;
- while (j > startndx && arr[j - 1] > arr[j]) {
- Swap(&arr[j],&arr[j-1]);
+ while (j > startndx && arr[j - 1] > arr[j]) {
+ Swap(&arr[j],&arr[j-1]);
pswap(&perm[j],&perm[j-1]);
- j--;
- }
+ j--;
+ }
- }
+ }
}
if(direction<0){
- for (i = startndx; i <=endndx; i++) {
- j = i;
+ for (i = startndx; i <=endndx; i++) {
+ j = i;
- while (j > startndx && arr[j - 1] < arr[j]) {
- Swap(&arr[j],&arr[j-1]);
+ while (j > startndx && arr[j - 1] < arr[j]) {
+ Swap(&arr[j],&arr[j-1]);
pswap(&perm[j],&perm[j-1]);
- j--;
- }
+ j--;
+ }
- }
+ }
}
}
int BinarySearch (real *array, int low, int high, real key,int direction)
{
-int mid, max, min;
-max=high+2;
-min=low+1;
-
-//Iterative implementation
-
-if (direction>=0){
-while (max-min>1){
- mid=(min+max)>>1;
- if(key<array[mid-1]) max=mid;
- else min=mid;
-}
-return min;
+ int mid, max, min;
+ max=high+2;
+ min=low+1;
+
+/*Iterative implementation*/
+
+ if (direction>=0){
+ while (max-min>1){
+ mid=(min+max)>>1;
+ if(key<array[mid-1]) max=mid;
+ else min=mid;
+ }
+ return min;
+ }
+
+ else if (direction<0){
+ while(max-min>1){
+ mid=(min+max)>>1;
+ if(key>array[mid-1]) max=mid;
+ else min=mid;
+ }
+ return min-1;
+
+ }/*end -ifelse direction*/
}
-else if (direction<0){
-while(max-min>1){
- mid=(min+max)>>1;
- if(key>array[mid-1]) max=mid;
- else min=mid;
-}
-return min-1;
-
- }//end -ifelse direction
-}
-
-int start_binsearch(real *array, int *perm, int low, int high, real key, int direction)
+int start_binsearch(real *array, int *perm, int low, int high,
+ real key, int direction)
{
insertionSort(array,perm,low,high,direction);
return BinarySearch(array,low,high,key,direction);
}
+
+int LinearSearch (double *array,int startindx, int stopindx,
+ double key,int *count, int direction)
+{
+ /*Iterative implementation - assume elements sorted*/
+ int i;
+ int keyindex;
+
+ if(direction>=0){
+ keyindex=startindx;
+ for (i=startindx;i<=stopindx;i++){
+ (*count)++;
+ if(array[i]>key) {
+ keyindex=i-1;
+ return keyindex;
+ }
+ }
+ }
+ else if(direction<0 ){
+ keyindex=stopindx;
+ for(i=stopindx;i>=startindx;i--){
+ (*count)++;
+ if (array[i]>key){
+ keyindex=i+1;
+ return keyindex;
+ }
+ }
+ }
+
+ else
+ fprintf(stderr,"Error: startindex=stopindex=%d\n",startindx);
+
+}
+
-#include <types/simple.h>
-void rangeArray(int *ar,int size);
-void insertionSort(real *ar, int *perm, int start, int end, int direction);
-int BinarySearch(real *ar, int start, int end, real key, int direction);
-int start_binsearch(real *array, int *perm, int low, int high, real key, int direction);
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ * $Id: densorder.c,v 0.9
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.0
+ *
+ * Copyright (c) 1991-2001
+ * BIOSON Research Institute, Dept. of Biophysical Chemistry
+ * University of Groningen, The Netherlands
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
+ *
+ * And Hey:
+ * Gyas ROwers Mature At Cryogenic Speed
+ */
+
+#ifndef _binsearch_h
+#define _binsearch_h
+
+#include "types/simple.h"
+
+extern void rangeArray(int *ar,int size);
+
+extern void insertionSort(real *ar, int *perm, int start, int end, int direction);
+
+extern int BinarySearch(real *ar, int start, int end, real key, int direction);
+
+extern int start_binsearch(real *array, int *perm, int low, int high,
+ real key, int direction);
+
+extern int LinearSearch(double *array,int startindx, int stopindx,
+ double key,int *count, int direction);
+
+#endif
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ * $Id: densorder.c,v 0.9
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.0
+ *
+ * Copyright (c) 1991-2001
+ * BIOSON Research Institute, Dept. of Biophysical Chemistry
+ * University of Groningen, The Netherlands
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
+ *
+ * And Hey:
+ * Gyas ROwers Mature At Cryogenic Speed
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
/* dens_filter.c
* Routines for Filters and convolutions
*/
#define EXP(x) (expf(x))
#endif
-// Modulo function assuming y>0 but with arbitrary INTEGER x
+/* Modulo function assuming y>0 but with arbitrary INTEGER x */
static int MOD(int x, int y){
-if (x<0) x=x+y;
-return (mod(x,y));
+ if (x<0) x=x+y;
+ return (mod(x,y));
}
int i, j, k;
real *out;
snew(out,dataSize);
- // check validity of params
+ /* check validity of params */
if(!x || !kernel) return FALSE;
if(dataSize <=0 || kernelSize <= 0) return FALSE;
- // start convolution from out[kernelSize-1] to out[dataSize-1] (last)
+ /* start convolution from out[kernelSize-1] to out[dataSize-1] (last) */
for(i = kernelSize-1; i < dataSize; ++i)
{
for(j = i, k = 0; k < kernelSize; --j, ++k)
out[i] += x[j] * kernel[k];
}
- // convolution from out[0] to out[kernelSize-2]
+ /* convolution from out[0] to out[kernelSize-2] */
for(i = 0; i < kernelSize - 1; ++i)
{
for(j = i, k = 0; j >= 0; --j, ++k)
return TRUE;
}
-//Assuming kernel is shorter than x
+/* Assuming kernel is shorter than x */
gmx_bool periodic_convolution(int datasize, real *x, int kernelsize, real *kernel)
{
- if (!x || !kernel) return FALSE;
- if (kernelsize<=0|| datasize<=0|| kernelsize > datasize) return FALSE;
+ if (!x || !kernel) return FALSE;
+ if (kernelsize<=0|| datasize<=0|| kernelsize > datasize) return FALSE;
int i,j,k,nj;
real *filtered;
for(i=0;(i<datasize); i++){
for(j=0; (j<kernelsize); j++){
- filtered[i] += kernel[j]*x[MOD((i-j),datasize)];
+ filtered[i] += kernel[j]*x[MOD((i-j),datasize)];
}
}
- for(i=0; i<datasize; i++) x[i]=filtered[i];
- sfree(filtered);
+ for(i=0; i<datasize; i++) x[i]=filtered[i];
+ sfree(filtered);
- return TRUE;
+ return TRUE;
}
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ * $Id: densorder.c,v 0.9
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.0
+ *
+ * Copyright (c) 1991-2001
+ * BIOSON Research Institute, Dept. of Biophysical Chemistry
+ * University of Groningen, The Netherlands
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
+ *
+ * And Hey:
+ * Gyas ROwers Mature At Cryogenic Speed
+ */
+
#ifndef _dens_filter_h
#define _dens_filter_h
#include "types/simple.h"
-gmx_bool convolution(int dataSize, real* in, int kernelSize, real* kernel);
-gmx_bool periodic_convolution(int dsize, real *in, int ksize, real* kernel);
-void gausskernel(real *out, int size, real var);
+extern gmx_bool convolution(int dataSize, real* in, int kernelSize,
+ real* kernel);
+extern gmx_bool periodic_convolution(int dsize, real *in, int ksize,
+ real* kernel);
+extern void gausskernel(real *out, int size, real var);
#endif
snew(dum,ma+1);
for(i=1; (i<ma+1); i++) {
lista[i] = i;
- ia[i] = 1; //fixed bug B.S.S 19/11
+ ia[i] = 1; /* fixed bug B.S.S 19/11 */
snew(covar[i],ma+1);
snew(alpha[i],ma+1);
}
-/*
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
* $Id: densorder.c,v 0.9
*
* This source code is part of
static void center_coords(t_atoms *atoms,matrix box,rvec x0[],int axis)
{
- int i,m;
- real tmass,mm;
- rvec com,shift,box_center;
+ int i,m;
+ real tmass,mm;
+ rvec com,shift,box_center;
- tmass = 0;
- clear_rvec(com);
- for(i=0; (i<atoms->nr); i++) {
- mm = atoms->atom[i].m;
- tmass += mm;
+ tmass = 0;
+ clear_rvec(com);
+ for(i=0; (i<atoms->nr); i++)
+ {
+ mm = atoms->atom[i].m;
+ tmass += mm;
+ for(m=0; (m<DIM); m++)
+ com[m] += mm*x0[i][m];
+ }
for(m=0; (m<DIM); m++)
- com[m] += mm*x0[i][m];
- }
- for(m=0; (m<DIM); m++)
- com[m] /= tmass;
- calc_box_center(ecenterDEF,box,box_center);
- rvec_sub(box_center,com,shift);
- shift[axis] -= box_center[axis];
+ com[m] /= tmass;
+ calc_box_center(ecenterDEF,box,box_center);
+ rvec_sub(box_center,com,shift);
+ shift[axis] -= box_center[axis];
- for(i=0; (i<atoms->nr); i++)
- rvec_dec(x0[i],shift);
+
for(i=0; (i<atoms->nr); i++)
+ rvec_dec(x0[i],shift);
}
* grpn - group number in index
*/
t_trxstatus *status;
- gmx_rmpbc_t gpbc=NULL;
+ gmx_rmpbc_t gpbc=NULL;
matrix box; /* Box - 3x3 -each step*/
rvec *x0; /* List of Coord without PBC*/
int natoms,i,j,k,n, /* loop indices, checks etc*/
- ax1=0,ax2=0, /* tangent directions */
- framenr=0, /* frame number in trajectory*/
- slicex, slicey, slicez; //slice # of x y z position
- real ***Densslice; // Density-slice in one frame
- real dscale; //physical scaling factor
- real t,x,y,z; // time and coordinates
- *tblock=0;/* blocknr in block average - initialise to 0*/
+ ax1=0,ax2=0, /* tangent directions */
+ framenr=0, /* frame number in trajectory*/
+ slicex, slicey, slicez; /*slice # of x y z position */
+ real ***Densslice; /* Density-slice in one frame*/
+ real dscale; /*physical scaling factor*/
+ real t,x,y,z; /* time and coordinates*/
rvec bbww;
+ *tblock=0;/* blocknr in block average - initialise to 0*/
/* Axis: X=0, Y=1,Z=2 */
switch(axis)
{
ax1=ZZ; ax2=XX; /* Surface : XZ*/
break;
case 2:
- ax1 = XX; ax2 = YY; /* Surface XY*/
- break;
+ ax1 = XX; ax2 = YY; /* Surface XY*/
+ break;
default:
- gmx_fatal(FARGS,"Invalid axes. Terminating\n");
+ gmx_fatal(FARGS,"Invalid axes. Terminating\n");
}
if( (natoms= read_first_x(oenv,&status,fn,&t,&x0,box)==0))
- gmx_fatal(FARGS, "Could not read coordinates from file"); // Open trajectory for read
+ gmx_fatal(FARGS, "Could not read coordinates from file"); /* Open trajectory for read*/
*zslices=1+FLOOR(box[axis][axis]/bwz);
*yslices=1+FLOOR(box[ax2][ax2]/bw);
*xslices=1+FLOOR(box[ax1][ax1]/bw);
- if(bps1d){
+ if(bps1d)
+ {
if(*xslices<*yslices) *xslices=1;
else *yslices=1;
- }
+ }
fprintf(stderr,
- "\nDividing the box in %5d x %5d x %5d slices with binw %f along axis %d\n",*xslices,*yslices,*zslices,bw,axis );
+ "\nDividing the box in %5d x %5d x %5d slices with binw %f along axis %d\n",*xslices,*yslices,*zslices,bw,axis );
/****Start trajectory processing***/
- //Initialize Densdevel and PBC-remove
+ /*Initialize Densdevel and PBC-remove*/
gpbc=gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box);
*Densdevel=NULL;
- do {
- bbww[XX] = box[ax1][ax1]/ *xslices;
- bbww[YY] = box[ax2][ax2]/ *yslices;
- bbww[ZZ] = box[axis][axis]/ *zslices;
- gmx_rmpbc(gpbc,top->atoms.nr,box,x0);
- //Reset Densslice every nsttblock steps
- if ( framenr % nsttblock==0 ){
+ do
+ {
+ bbww[XX] = box[ax1][ax1]/ *xslices;
+ bbww[YY] = box[ax2][ax2]/ *yslices;
+ bbww[ZZ] = box[axis][axis]/ *zslices;
+ gmx_rmpbc(gpbc,top->atoms.nr,box,x0);
+ /*Reset Densslice every nsttblock steps*/
+ if ( framenr % nsttblock==0 )
+ {
snew(Densslice,*xslices);
- for (i=0;i<*xslices;i++) {
- snew(Densslice[i],*yslices);
- for (j=0;j<*yslices;j++){
- snew(Densslice[i][j],*zslices);
- }
- }
+ for (i=0;i<*xslices;i++)
+ {
+ snew(Densslice[i],*yslices);
+ for (j=0;j<*yslices;j++)
+ {
+ snew(Densslice[i][j],*zslices);
+ }
+ }
- /*Allocate Memory to extra frame in Densdevel - rather stupid approach: *A single frame each time, although only every nsttblock steps.*/
+ /*Allocate Memory to extra frame in Densdevel - rather stupid approach: *A single frame each time, although only every nsttblock steps.*/
srenew(*Densdevel,*tblock+1);
}
center_coords(&top->atoms,box,x0,axis);
- for (j=0;j<gnx[0];j++) { /*Loop over all atoms in selected index*/
+ for (j=0;j<gnx[0];j++)
+ { /*Loop over all atoms in selected index*/
x=x0[index[0][j]][ax1];
y=x0[index[0][j]][ax2];
z=x0[index[0][j]][axis];
framenr++;
- if(framenr % nsttblock == 0){
- //Implicit incrementation of Densdevel via renewal of Densslice
- //only every nsttblock steps
+ if(framenr % nsttblock == 0)
+ {
+ /*Implicit incrementation of Densdevel via renewal of Densslice*/
+ /*only every nsttblock steps*/
(*Densdevel)[*tblock]=Densslice;
(*tblock)++;
}
- } while(read_next_x(oenv,status,&t,natoms,x0,box));
+ } while(read_next_x(oenv,status,&t,natoms,x0,box));
- //Free memory we no longer need and exit.
+ /*Free memory we no longer need and exit.*/
gmx_rmpbc_done(gpbc);
close_trj(status);
if (0)
{
- FILE *fp;
- fp = fopen("koko.xvg","w");
- for(j=0; (j<*zslices); j++) {
- fprintf(fp,"%5d",j);
- for(i=0; (i<*tblock); i++) {
- fprintf(fp," %10g",(*Densdevel)[i][9][1][j]);
- }
- fprintf(fp,"\n");
- }
- fclose(fp);
+ FILE *fp;
+ fp = fopen("koko.xvg","w");
+ for(j=0; (j<*zslices); j++)
+ {
+ fprintf(fp,"%5d",j);
+ for(i=0; (i<*tblock); i++)
+ {
+ fprintf(fp," %10g",(*Densdevel)[i][9][1][j]);
+ }
+ fprintf(fp,"\n");
+ }
+ fclose(fp);
}
}
static void printdensavg(char *fldfn, real ****Densmap, int xslices, int yslices, int zslices, int tdim)
{
/*Debug-filename and filehandle*/
- FILE *fldH;
- int n,i,j,k;
- real totdens=0;
- fldH=ffopen(fldfn,"w");
- fprintf(fldH,"%i %i %i %i\n",tdim, xslices,yslices,zslices);
- for(n=0;n<tdim;n++){
- for(i=0;i<xslices;i++){
- for (j=0;j<yslices;j++){
- for (k=0;k<zslices;k++){
- fprintf(fldH,"%i %i %i %f\n",i,j,k,Densmap[n][i][j][k]);
- totdens+=(Densmap[n][i][j][k]);
- }
- }
- }
- }
-totdens/=(xslices*yslices*zslices*tdim);
-fprintf(stderr,"Total density [kg/m^3] %8f",totdens);
-ffclose(fldH);
+ FILE *fldH;
+ int n,i,j,k;
+ real totdens=0;
+ fldH=ffopen(fldfn,"w");
+ fprintf(fldH,"%i %i %i %i\n",tdim, xslices,yslices,zslices);
+ for(n=0;n<tdim;n++)
+ {
+ for(i=0;i<xslices;i++)
+ {
+ for (j=0;j<yslices;j++)
+ {
+ for (k=0;k<zslices;k++)
+ {
+ fprintf(fldH,"%i %i %i %f\n",i,j,k,Densmap[n][i][j][k]);
+ totdens+=(Densmap[n][i][j][k]);
+ }
+ }
+ }
+ }
+ totdens/=(xslices*yslices*zslices*tdim);
+ fprintf(stderr,"Total density [kg/m^3] %8f",totdens);
+ ffclose(fldH);
}
-static void outputfield(char *fldfn, real ****Densmap, int xslices, int yslices, int zslices, int tdim)
+static void outputfield(const char *fldfn, real ****Densmap,
+ int xslices, int yslices, int zslices, int tdim)
{
/*Debug-filename and filehandle*/
- FILE *fldH;
- int n,i,j,k;
- int dim[4]={tdim,xslices,yslices,zslices};
- real totdens=0;
- fldH=ffopen(fldfn,"w");
- fwrite(dim,sizeof(int),4,fldH);
- for(n=0;n<tdim;n++){
- for(i=0;i<xslices;i++){
- for (j=0;j<yslices;j++){
- for (k=0;k<zslices;k++){
- fwrite(&(Densmap[n][i][j][k]),sizeof(real),1,fldH);
- totdens+=(Densmap[n][i][j][k]);
- }
- }
- }
- }
-totdens/=(xslices*yslices*zslices*tdim);
-fprintf(stderr,"Total density [kg/m^3] %8f",totdens);
-ffclose(fldH);
+ FILE *fldH;
+ int n,i,j,k;
+ int dim[4];
+ real totdens=0;
+
+ dim[0] = tdim;
+ dim[1] = xslices;
+ dim[2] = yslices;
+ dim[3] = zslices;
+
+ fldH=ffopen(fldfn,"w");
+ fwrite(dim,sizeof(int),4,fldH);
+ for(n=0;n<tdim;n++)
+ {
+ for(i=0;i<xslices;i++)
+ {
+ for (j=0;j<yslices;j++)
+ {
+ for (k=0;k<zslices;k++)
+ {
+ fwrite(&(Densmap[n][i][j][k]),sizeof(real),1,fldH);
+ totdens+=(Densmap[n][i][j][k]);
+ }
+ }
+ }
+ }
+ totdens/=(xslices*yslices*zslices*tdim);
+ fprintf(stderr,"Total density [kg/m^3] %8f",totdens);
+ ffclose(fldH);
}
static void periodic_running_average(int npoints,real *x,int naver)
{
- int i,j,nj;
- real *aver;
+ int i,j,nj;
+ real *aver;
- snew(aver,npoints);
- for(i=0; (i<npoints); i++) {
- nj = 0;
- for(j=i-naver/2; (j<i+naver/2); j++) {
- aver[i] += x[(j+npoints) % npoints];
- nj++;
+ snew(aver,npoints);
+ for(i=0; (i<npoints); i++)
+ {
+ nj = 0;
+ for(j=i-naver/2; (j<i+naver/2); j++)
+ {
+ aver[i] += x[(j+npoints) % npoints];
+ nj++;
+ }
+ aver[i] /= nj;
}
- aver[i] /= nj;
- }
- for(i=0; (i<npoints); i++)
- x[i] = aver[i];
- sfree(aver);
+ for(i=0; (i<npoints); i++)
+ x[i] = aver[i];
+ sfree(aver);
}
static void filterdensmap(real ****Densmap, int xslices, int yslices, int zslices,int tblocks,int ftsize)
{
- real *kernel;
- real *output;
- real std,var;
- int i,j,k,n, order;
- order=ftsize/2;
- std=((real)order/2.0);
- var=std*std;
- snew(kernel,ftsize);
- gausskernel(kernel,ftsize,var);
- for(n=0;n<tblocks;n++){
- for(i=0;i<xslices;i++){
- for (j=0;j<yslices;j++){
- periodic_convolution(zslices,Densmap[n][i][j],ftsize,kernel);
- if (0) periodic_running_average(zslices,Densmap[n][i][j],ftsize);
- }
- }
- }
+ real *kernel;
+ real *output;
+ real std,var;
+ int i,j,k,n, order;
+ order=ftsize/2;
+ std=((real)order/2.0);
+ var=std*std;
+ snew(kernel,ftsize);
+ gausskernel(kernel,ftsize,var);
+ for(n=0;n<tblocks;n++)
+ {
+ for(i=0;i<xslices;i++)
+ {
+ for (j=0;j<yslices;j++)
+ {
+ periodic_convolution(zslices,Densmap[n][i][j],ftsize,kernel);
+ }
+ }
+ }
}
static void interfaces_txy (real ****Densmap, int xslices, int yslices, int zslices,
- int tblocks, real binwidth,int method, real dens1, real dens2, t_interf ****intf1, t_interf ****intf2,const output_env_t oenv)
+ int tblocks, real binwidth,int method,
+ real dens1, real dens2, t_interf ****intf1,
+ t_interf ****intf2,const output_env_t oenv)
{
-
-//Returns two pointers to 3D arrays of t_interf structs containing (position,thickness) of the interface(s)
+ /*Returns two pointers to 3D arrays of t_interf structs containing (position,thickness) of the interface(s)*/
FILE *xvg;
- real *zDensavg; // zDensavg[z]
+ real *zDensavg; /* zDensavg[z]*/
int i,j,k,n;
int xysize;
int ndx1, ndx2, deltandx, *zperm;
real densmid, densl, densr, alpha, pos, spread;
real splitpoint, startpoint, endpoint;
real *sigma1, *sigma2;
+ real beginfit1[4];
+ real beginfit2[4];
+ real *fit1=NULL,*fit2=NULL;
+ const real *avgfit1;
+ const real *avgfit2;
const real onehalf= 1.00/2.00;
- t_interf ***int1=NULL,***int2=NULL; //Interface matrices [t][x,y] - last index in row-major order
- /*Create int1(t,xy) and int2(t,xy) arrays with correct number of interf_t elements*/
+ t_interf ***int1=NULL,***int2=NULL; /*Interface matrices [t][x,y] - last index in row-major order*/
+ /*Create int1(t,xy) and int2(t,xy) arrays with correct number of interf_t elements*/
xysize=xslices*yslices;
snew(int1,tblocks);
snew(int2,tblocks);
- for (i=0;i<tblocks;i++){
+ for (i=0;i<tblocks;i++)
+ {
snew(int1[i],xysize);
snew(int2[i],xysize);
- for(j=0;j<xysize;j++){
+ for(j=0;j<xysize;j++)
+ {
snew(int1[i][j],1);
snew(int2[i][j],1);
init_interf(int1[i][j]);
}
}
-if(method==methBISECT){
- densmid= onehalf*(dens1+dens2);
- snew(zperm,zslices);
- for(n=0;n<tblocks;n++){
- for(i=0;i<xslices;i++){
- for (j=0;j<yslices;j++){
- rangeArray(zperm,zslices); //reset permutation array to identity
- /*Binsearch returns slice-nr where the order param is <= setpoint sgmid*/
- ndx1=start_binsearch(Densmap[n][i][j],zperm,0,zslices/2-1,densmid,1);
- ndx2=start_binsearch(Densmap[n][i][j],zperm,zslices/2,zslices-1,densmid,-1);
+ if(method==methBISECT)
+ {
+ densmid= onehalf*(dens1+dens2);
+ snew(zperm,zslices);
+ for(n=0;n<tblocks;n++)
+ {
+ for(i=0;i<xslices;i++)
+ {
+ for (j=0;j<yslices;j++)
+ {
+ rangeArray(zperm,zslices); /*reset permutation array to identity*/
+ /*Binsearch returns slice-nr where the order param is <= setpoint sgmid*/
+ ndx1=start_binsearch(Densmap[n][i][j],zperm,0,zslices/2-1,densmid,1);
+ ndx2=start_binsearch(Densmap[n][i][j],zperm,zslices/2,zslices-1,densmid,-1);
- /* Linear interpolation (for use later if time allows)
- * rho_1s= Densmap[n][i][j][zperm[ndx1]]
- * rho_1e =Densmap[n][i][j][zperm[ndx1+1]] - in worst case might be far off
- * rho_2s =Densmap[n][i][j][zperm[ndx2+1]]
- * rho_2e =Densmap[n][i][j][zperm[ndx2]]
- * For 1st interface we have:
- densl= Densmap[n][i][j][zperm[ndx1]];
- densr= Densmap[n][i][j][zperm[ndx1+1]];
- alpha=(densmid-densl)/(densr-densl);
- deltandx=zperm[ndx1+1]-zperm[ndx1];
-
- if(debug){
- printf("Alpha, Deltandx %f %i\n", alpha,deltandx);
- }
- if(abs(alpha)>1.0 || abs(deltandx)>3){
- pos=zperm[ndx1];
- spread=-1;
- }
- else {
- pos=zperm[ndx1]+alpha*deltandx;
- spread=binwidth*deltandx;
- }
- * For the 2nd interface can use the same formulation, since alpha should become negative ie:
- * alpha=(densmid-Densmap[n][i][j][zperm[ndx2]])/(Densmap[n][i][j][zperm[nxd2+1]]-Densmap[n][i][j][zperm[ndx2]]);
- * deltandx=zperm[ndx2+1]-zperm[ndx2];
- * pos=zperm[ndx2]+alpha*deltandx; */
-
- //After filtering we use the direct approach
- int1[n][j+(i*yslices)]->Z=(zperm[ndx1]+onehalf)*binwidth;
- int1[n][j+(i*yslices)]->t=binwidth;
- int2[n][j+(i*yslices)]->Z=(zperm[ndx2]+onehalf)*binwidth;
- int2[n][j+(i*yslices)]->t=binwidth;
- }
- }
- }
-}
-
-if(method==methFUNCFIT){
- //Assume a box divided in 2 along midpoint of z for starters
- startpoint=0.0;
- endpoint = binwidth*zslices;
- splitpoint = (startpoint+endpoint)/2.0;
- //Initial fit proposals
- real beginfit1[4] = {dens1, dens2, (splitpoint/2), 0.5};
- real beginfit2[4] = {dens2, dens1, (3*splitpoint/2), 0.5};
- real *fit1=NULL,*fit2=NULL;
-
- snew(zDensavg,zslices);
- snew(sigma1,zslices);
- snew(sigma2,zslices);
-
- for(k=0;k<zslices;k++) sigma1[k]=sigma2[k]=1;
- //Calculate average density along z - avoid smoothing by using coarse-grained-mesh
- for(k=0;k<zslices;k++){
- for(n=0;n<tblocks;n++){
- for (i=0;i<xslices;i++){
- for(j=0;j<yslices;j++){
- zDensavg[k]+=(Densmap[n][i][j][k]/(xslices*yslices*tblocks));
- }
- }
- }
- }
-
- if(debug){
- xvg=xvgropen("DensprofileonZ.xvg", "Averaged Densityprofile on Z","z[nm]","Density[kg/m^3]",oenv);
- for(k=0;k<zslices;k++) fprintf(xvg, "%4f.3 %8f.4\n", k*binwidth,zDensavg[k]);
- fclose(xvg);
- }
+ /* Linear interpolation (for use later if time allows)
+ * rho_1s= Densmap[n][i][j][zperm[ndx1]]
+ * rho_1e =Densmap[n][i][j][zperm[ndx1+1]] - in worst case might be far off
+ * rho_2s =Densmap[n][i][j][zperm[ndx2+1]]
+ * rho_2e =Densmap[n][i][j][zperm[ndx2]]
+ * For 1st interface we have:
+ densl= Densmap[n][i][j][zperm[ndx1]];
+ densr= Densmap[n][i][j][zperm[ndx1+1]];
+ alpha=(densmid-densl)/(densr-densl);
+ deltandx=zperm[ndx1+1]-zperm[ndx1];
+
+ if(debug){
+ printf("Alpha, Deltandx %f %i\n", alpha,deltandx);
+ }
+ if(abs(alpha)>1.0 || abs(deltandx)>3){
+ pos=zperm[ndx1];
+ spread=-1;
+ }
+ else {
+ pos=zperm[ndx1]+alpha*deltandx;
+ spread=binwidth*deltandx;
+ }
+ * For the 2nd interface can use the same formulation, since alpha should become negative ie:
+ * alpha=(densmid-Densmap[n][i][j][zperm[ndx2]])/(Densmap[n][i][j][zperm[nxd2+1]]-Densmap[n][i][j][zperm[ndx2]]);
+ * deltandx=zperm[ndx2+1]-zperm[ndx2];
+ * pos=zperm[ndx2]+alpha*deltandx; */
+
+ /*After filtering we use the direct approach */
+ int1[n][j+(i*yslices)]->Z=(zperm[ndx1]+onehalf)*binwidth;
+ int1[n][j+(i*yslices)]->t=binwidth;
+ int2[n][j+(i*yslices)]->Z=(zperm[ndx2]+onehalf)*binwidth;
+ int2[n][j+(i*yslices)]->t=binwidth;
+ }
+ }
+ }
+ }
+
+ if(method==methFUNCFIT)
+ {
+ /*Assume a box divided in 2 along midpoint of z for starters*/
+ startpoint=0.0;
+ endpoint = binwidth*zslices;
+ splitpoint = (startpoint+endpoint)/2.0;
+ /*Initial fit proposals*/
+ beginfit1[0] = dens1;
+ beginfit1[1] = dens2;
+ beginfit1[2] = (splitpoint/2);
+ beginfit1[3] = 0.5;
+
+ beginfit2[0] = dens2;
+ beginfit2[1] = dens1;
+ beginfit2[2] = (3*splitpoint/2);
+ beginfit2[3] = 0.5;
+
+ snew(zDensavg,zslices);
+ snew(sigma1,zslices);
+ snew(sigma2,zslices);
+
+ for(k=0;k<zslices;k++) sigma1[k]=sigma2[k]=1;
+ /*Calculate average density along z - avoid smoothing by using coarse-grained-mesh*/
+ for(k=0;k<zslices;k++)
+ {
+ for(n=0;n<tblocks;n++)
+ {
+ for (i=0;i<xslices;i++)
+ {
+ for(j=0;j<yslices;j++)
+ {
+ zDensavg[k]+=(Densmap[n][i][j][k]/(xslices*yslices*tblocks));
+ }
+ }
+ }
+ }
+
+ if(debug){
+ xvg=xvgropen("DensprofileonZ.xvg", "Averaged Densityprofile on Z","z[nm]","Density[kg/m^3]",oenv);
+ for(k=0;k<zslices;k++) fprintf(xvg, "%4f.3 %8f.4\n", k*binwidth,zDensavg[k]);
+ fclose(xvg);
+ }
- /*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/
+ /*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/
- //Fit 1st half of box
- do_lmfit(zslices, zDensavg,sigma1,binwidth,NULL,startpoint,splitpoint,oenv,FALSE,effnERF,beginfit1,3);
- //Fit 2nd half of box
- do_lmfit(zslices ,zDensavg,sigma2,binwidth,NULL,splitpoint,endpoint,oenv,FALSE,effnERF,beginfit2,3);
+ /*Fit 1st half of box*/
+ do_lmfit(zslices, zDensavg,sigma1,binwidth,NULL,startpoint,splitpoint,oenv,FALSE,effnERF,beginfit1,3);
+ /*Fit 2nd half of box*/
+ do_lmfit(zslices ,zDensavg,sigma2,binwidth,NULL,splitpoint,endpoint,oenv,FALSE,effnERF,beginfit2,3);
- //Initialise the const arrays for storing the average fit parameters
- const real * const avgfit1=beginfit1;
- const real * const avgfit2=beginfit2;
+ /*Initialise the const arrays for storing the average fit parameters*/
+ avgfit1=beginfit1;
+ avgfit2=beginfit2;
- //Now do fit over each x y and t slice to get Zint(x,y,t) - loop is very large, we potentially should average over time directly
- for(n=0;n<tblocks;n++){
- for(i=0;i<xslices;i++){
- for (j=0;j<yslices;j++){
- //Reinitialise fit for each mesh-point
- srenew(fit1,4);
- srenew(fit2,4);
- for (k=0;k<4;k++){
- fit1[k]=avgfit1[k];
- fit2[k]=avgfit2[k];
- }
- //Now fit and store in structures in row-major order int[n][i][j]
- do_lmfit(zslices,Densmap[n][i][j],sigma1,binwidth,NULL,startpoint,splitpoint,oenv,FALSE,effnERF,fit1,1);
- int1[n][j+(yslices*i)]->Z=fit1[2];
- int1[n][j+(yslices*i)]->t=fit1[3];
- do_lmfit(zslices,Densmap[n][i][j],sigma2,binwidth,NULL,splitpoint,endpoint,oenv,FALSE,effnERF,fit2,2);
- int2[n][j+(yslices*i)]->Z=fit2[2];
- int2[n][j+(yslices*i)]->t=fit2[3];
- }
- }
- }
-}
+ /*Now do fit over each x y and t slice to get Zint(x,y,t) - loop is very large, we potentially should average over time directly*/
+ for(n=0;n<tblocks;n++)
+ {
+ for(i=0;i<xslices;i++)
+ {
+ for (j=0;j<yslices;j++)
+ {
+ /*Reinitialise fit for each mesh-point*/
+ srenew(fit1,4);
+ srenew(fit2,4);
+ for (k=0;k<4;k++)
+ {
+ fit1[k]=avgfit1[k];
+ fit2[k]=avgfit2[k];
+ }
+ /*Now fit and store in structures in row-major order int[n][i][j]*/
+ do_lmfit(zslices,Densmap[n][i][j],sigma1,binwidth,NULL,startpoint,splitpoint,oenv,FALSE,effnERF,fit1,1);
+ int1[n][j+(yslices*i)]->Z=fit1[2];
+ int1[n][j+(yslices*i)]->t=fit1[3];
+ do_lmfit(zslices,Densmap[n][i][j],sigma2,binwidth,NULL,splitpoint,endpoint,oenv,FALSE,effnERF,fit2,2);
+ int2[n][j+(yslices*i)]->Z=fit2[2];
+ int2[n][j+(yslices*i)]->t=fit2[3];
+ }
+ }
+ }
+ }
-*intf1=int1;
-*intf2=int2;
+ *intf1=int1;
+ *intf2=int2;
}
static void writesurftoxpms(t_interf ***surf1,t_interf ***surf2, int tblocks,int xbins, int ybins, int zbins, real bw,real bwz, char **outfiles,int maplevels )
{
- char numbuf[12];
- int n, i, j;
- real **profile1, **profile2;
- real max1, max2, min1, min2, *xticks, *yticks;
- t_rgb lo={0,0,0};
- t_rgb hi={1,1,1};
- FILE *xpmfile1, *xpmfile2;
+ char numbuf[12];
+ int n, i, j;
+ real **profile1, **profile2;
+ real max1, max2, min1, min2, *xticks, *yticks;
+ t_rgb lo={0,0,0};
+ t_rgb hi={1,1,1};
+ FILE *xpmfile1, *xpmfile2;
/*Prepare xpm structures for output*/
-//Allocate memory to tick's and matrices
-snew (xticks,xbins+1);
-snew (yticks,ybins+1);
+/*Allocate memory to tick's and matrices*/
+ snew (xticks,xbins+1);
+ snew (yticks,ybins+1);
-profile1=mk_matrix(xbins,ybins,FALSE);
-profile2=mk_matrix(xbins,ybins,FALSE);
+ profile1=mk_matrix(xbins,ybins,FALSE);
+ profile2=mk_matrix(xbins,ybins,FALSE);
-for (i=0;i<xbins+1;i++) xticks[i]+=bw;
-for (j=0;j<ybins+1;j++) yticks[j]+=bw;
-
-xpmfile1 = ffopen(outfiles[0],"w");
-xpmfile2 = ffopen(outfiles[1],"w");
-
-max1=max2=0.0;
-min1=min2=zbins*bwz;
-
-for(n=0;n<tblocks;n++){
-sprintf(numbuf,"tblock: %4i",n);
-//Filling matrices for inclusion in xpm-files
- for(i=0;i<xbins;i++){
- for(j=0;j<ybins;j++){
+ for (i=0;i<xbins+1;i++) xticks[i]+=bw;
+ for (j=0;j<ybins+1;j++) yticks[j]+=bw;
+
+ xpmfile1 = ffopen(outfiles[0],"w");
+ xpmfile2 = ffopen(outfiles[1],"w");
+
+ max1=max2=0.0;
+ min1=min2=zbins*bwz;
+
+ for(n=0;n<tblocks;n++)
+ {
+ sprintf(numbuf,"tblock: %4i",n);
+/*Filling matrices for inclusion in xpm-files*/
+ for(i=0;i<xbins;i++)
+ {
+ for(j=0;j<ybins;j++)
+ {
profile1[i][j]=(surf1[n][j+ybins*i])->Z;
profile2[i][j]=(surf2[n][j+ybins*i])->Z;
- //Finding max and min values
+ /*Finding max and min values*/
if(profile1[i][j]>max1) max1=profile1[i][j];
if(profile1[i][j]<min1) min1=profile1[i][j];
if(profile2[i][j]>max2) max2=profile2[i][j];
if(profile2[i][j]<min2) min2=profile2[i][j];
- }
- }
+ }
+ }
- write_xpm(xpmfile1,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile1,min1,max1,lo,hi,&maplevels);
- write_xpm(xpmfile2,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile2,min2,max2,lo,hi,&maplevels);
-}
+ write_xpm(xpmfile1,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile1,min1,max1,lo,hi,&maplevels);
+ write_xpm(xpmfile2,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile2,min2,max2,lo,hi,&maplevels);
+ }
-ffclose(xpmfile1);
-ffclose(xpmfile2);
+ ffclose(xpmfile1);
+ ffclose(xpmfile2);
-sfree(profile1);
-sfree(profile2);
-sfree(xticks);
-sfree(yticks);
+ sfree(profile1);
+ sfree(profile2);
+ sfree(xticks);
+ sfree(yticks);
}
-static void writeraw(t_interf ***int1, t_interf ***int2, int tblocks,int xbins, int ybins,char **fnms){
+static void writeraw(t_interf ***int1, t_interf ***int2, int tblocks,int xbins, int ybins,char **fnms)
+{
FILE *raw1, *raw2;
- int i,j,n;
+ int i,j,n;
raw1=ffopen(fnms[0],"w");
raw2=ffopen(fnms[1],"w");
fprintf(raw2,"#Legend: nt nx ny\n#Xbin Ybin Z t\n");
fprintf(raw1,"%i %i %i\n", tblocks,xbins,ybins);
fprintf(raw2,"%i %i %i\n", tblocks,xbins,ybins);
- for (n=0;n<tblocks;n++){
- for(i=0;i<xbins;i++){
- for(j=0;j<ybins;j++){
+ for (n=0;n<tblocks;n++)
+ {
+ for(i=0;i<xbins;i++)
+ {
+ for(j=0;j<ybins;j++)
+ {
fprintf(raw1,"%i %i %8.5f %6.4f\n",i,j,(int1[n][j+ybins*i])->Z,(int1[n][j+ybins*i])->t);
fprintf(raw2,"%i %i %8.5f %6.4f\n",i,j,(int2[n][j+ybins*i])->Z,(int2[n][j+ybins*i])->t);
}
int gmx_densorder(int argc,char *argv[])
{
- static const char *desc[] = {
- "A small program to reduce a two-phase density distribution",
- "along an axis, computed over a MD trajectory",
- "to 2D surfaces fluctuating in time, by a fit to",
- "a functional profile for interfacial densities",
- "A time-averaged spatial representation of the" ,
- "interfaces can be output with the option -tavg"
- };
+ static const char *desc[] = {
+ "A small program to reduce a two-phase density distribution",
+ "along an axis, computed over a MD trajectory",
+ "to 2D surfaces fluctuating in time, by a fit to",
+ "a functional profile for interfacial densities",
+ "A time-averaged spatial representation of the" ,
+ "interfaces can be output with the option -tavg"
+ };
- /* Extra arguments - but note how you always get the begin/end
- * options when running the program, without mentioning them here!
- */
+ /* Extra arguments - but note how you always get the begin/end
+ * options when running the program, without mentioning them here!
+ */
- output_env_t oenv;
- t_topology *top;
- char title[STRLEN],**grpname;
- int ePBC, *ngx;
- static real binw=0.2;
- static real binwz=0.05;
- static real dens1=0.00;
- static real dens2=1000.00;
- static int ftorder=0;
- static int nsttblock=100;
- static int axis= 2;
- static char *axtitle="Z";
- static int ngrps=1;
- atom_id **index; // Index list for single group
- int xslices, yslices, zslices, tblock;
- static gmx_bool bGraph=FALSE;
- static gmx_bool bCenter = FALSE;
- static gmx_bool bFourier=FALSE;
- static gmx_bool bRawOut=FALSE;
- static gmx_bool bOut=FALSE;
- static gmx_bool b1d=FALSE;
- static int nlevels=100;
-//Densitymap - Densmap[t][x][y][z]
- real ****Densmap=NULL;
-// Surfaces surf[t][surf_x,surf_y]
- t_interf ***surf1, ***surf2;
-
- static const char *meth[]={NULL,"bisect","functional",NULL};
- int eMeth;
-
- char **graphfiles, **rawfiles, **spectra; /* Filenames for xpm-surface maps, rawdata and powerspectra */
- int nfxpm,nfraw, nfspect; /* # files for interface maps and spectra = # interfaces */
+ output_env_t oenv;
+ t_topology *top;
+ char title[STRLEN],**grpname;
+ int ePBC, *ngx;
+ static real binw=0.2;
+ static real binwz=0.05;
+ static real dens1=0.00;
+ static real dens2=1000.00;
+ static int ftorder=0;
+ static int nsttblock=100;
+ static int axis= 2;
+ static char *axtitle="Z";
+ static int ngrps=1;
+ atom_id **index; /* Index list for single group*/
+ int xslices, yslices, zslices, tblock;
+ static gmx_bool bGraph=FALSE;
+ static gmx_bool bCenter = FALSE;
+ static gmx_bool bFourier=FALSE;
+ static gmx_bool bRawOut=FALSE;
+ static gmx_bool bOut=FALSE;
+ static gmx_bool b1d=FALSE;
+ static int nlevels=100;
+ /*Densitymap - Densmap[t][x][y][z]*/
+ real ****Densmap=NULL;
+ /* Surfaces surf[t][surf_x,surf_y]*/
+ t_interf ***surf1, ***surf2;
+
+ static const char *meth[]={NULL,"bisect","functional",NULL};
+ int eMeth;
+
+ char **graphfiles, **rawfiles, **spectra; /* Filenames for xpm-surface maps, rawdata and powerspectra */
+ int nfxpm,nfraw, nfspect; /* # files for interface maps and spectra = # interfaces */
-t_pargs pa[] = {
- { "-1d", FALSE, etBOOL, {&b1d},
- "Pseudo-1d interface geometry"},
- { "-bw",FALSE,etREAL,{&binw},
- "Binwidth of density distribution tangential to interface"},
- { "-bwn", FALSE, etREAL, {&binwz},
- "Binwidth of density distribution normal to interface"},
- { "-order", FALSE, etINT,{&ftorder},
- "Order of Gaussian filter, order 0 equates to NO filtering"},
- {"-axis", FALSE, etSTR,{&axtitle},
- "Axis Direction - X, Y or Z"},
- {"-method",FALSE ,etENUM,{meth},
- "Interface location method"},
- {"-d1", FALSE, etREAL,{&dens1},
- "Bulk density phase 1 (at small z)"},
- {"-d2", FALSE, etREAL,{&dens2},
- "Bulk density phase 2 (at large z)"},
- { "-tblock",FALSE,etINT,{&nsttblock},
- "Number of frames in one time-block average"},
- { "-nlevel", FALSE,etINT, {&nlevels},
- "Number of Height levels in 2D - XPixMaps"}
- };
-
-
- t_filenm fnm[] = {
- { efTPX, "-s", NULL, ffREAD }, /* this is for the topology */
- { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */
- { efNDX, "-n", NULL, ffREAD}, /* this is to select groups */
- { efDAT, "-o", "Density4D",ffOPTWR}, /* This is for outputting the entire 4D densityfield in binary format */
- { efOUT, "-or", NULL,ffOPTWRMULT}, /* This is for writing out the entire information in the t_interf arrays */
- { efXPM, "-og" ,"interface",ffOPTWRMULT},/* This is for writing out the interface meshes - one xpm-file per tblock*/
- { efOUT, "-Spect","intfspect",ffOPTWRMULT}, /* This is for the trajectory averaged Fourier-spectra*/
- };
+ t_pargs pa[] = {
+ { "-1d", FALSE, etBOOL, {&b1d},
+ "Pseudo-1d interface geometry"},
+ { "-bw",FALSE,etREAL,{&binw},
+ "Binwidth of density distribution tangential to interface"},
+ { "-bwn", FALSE, etREAL, {&binwz},
+ "Binwidth of density distribution normal to interface"},
+ { "-order", FALSE, etINT,{&ftorder},
+ "Order of Gaussian filter, order 0 equates to NO filtering"},
+ {"-axis", FALSE, etSTR,{&axtitle},
+ "Axis Direction - X, Y or Z"},
+ {"-method",FALSE ,etENUM,{meth},
+ "Interface location method"},
+ {"-d1", FALSE, etREAL,{&dens1},
+ "Bulk density phase 1 (at small z)"},
+ {"-d2", FALSE, etREAL,{&dens2},
+ "Bulk density phase 2 (at large z)"},
+ { "-tblock",FALSE,etINT,{&nsttblock},
+ "Number of frames in one time-block average"},
+ { "-nlevel", FALSE,etINT, {&nlevels},
+ "Number of Height levels in 2D - XPixMaps"}
+ };
+
+
+ t_filenm fnm[] = {
+ { efTPX, "-s", NULL, ffREAD }, /* this is for the topology */
+ { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */
+ { efNDX, "-n", NULL, ffREAD}, /* this is to select groups */
+ { efDAT, "-o", "Density4D",ffOPTWR}, /* This is for outputting the entire 4D densityfield in binary format */
+ { efOUT, "-or", NULL,ffOPTWRMULT}, /* This is for writing out the entire information in the t_interf arrays */
+ { efXPM, "-og" ,"interface",ffOPTWRMULT},/* This is for writing out the interface meshes - one xpm-file per tblock*/
+ { efOUT, "-Spect","intfspect",ffOPTWRMULT}, /* This is for the trajectory averaged Fourier-spectra*/
+ };
#define NFILE asize(fnm)
- CopyRight(stderr,argv[0]);
+ CopyRight(stderr,argv[0]);
- /* This is the routine responsible for adding default options,
- * calling the X/motif interface, etc. */
- parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,
- NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
+ /* This is the routine responsible for adding default options,
+ * calling the X/motif interface, etc. */
+ parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,
+ NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
-eMeth=nenum(meth);
-bFourier=opt2bSet("-Spect",NFILE,fnm);
-bRawOut=opt2bSet("-or",NFILE,fnm);
-bGraph=opt2bSet("-og",NFILE,fnm);
-bOut=opt2bSet("-o",NFILE,fnm);
-top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
-snew(grpname,ngrps);
-snew(index,ngrps);
-snew(ngx,ngrps);
+ eMeth=nenum(meth);
+ bFourier=opt2bSet("-Spect",NFILE,fnm);
+ bRawOut=opt2bSet("-or",NFILE,fnm);
+ bGraph=opt2bSet("-og",NFILE,fnm);
+ bOut=opt2bSet("-o",NFILE,fnm);
+ top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
+ snew(grpname,ngrps);
+ snew(index,ngrps);
+ snew(ngx,ngrps);
/* Calculate axis */
- axis = toupper(axtitle[0]) - 'X';
+ axis = toupper(axtitle[0]) - 'X';
-get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),ngrps,ngx,index,grpname);
+ get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),ngrps,ngx,index,grpname);
-density_in_time(ftp2fn(efTRX,NFILE,fnm),index,ngx,ngrps,binw,binwz,nsttblock,&Densmap,&xslices,&yslices,&zslices,&tblock,top,ePBC,axis,bCenter,b1d,oenv);
+ density_in_time(ftp2fn(efTRX,NFILE,fnm),index,ngx,ngrps,binw,binwz,nsttblock,&Densmap,&xslices,&yslices,&zslices,&tblock,top,ePBC,axis,bCenter,b1d,oenv);
-if(ftorder>0){
-filterdensmap(Densmap,xslices,yslices,zslices,tblock,2*ftorder+1);
-}
+ if(ftorder>0)
+ {
+ filterdensmap(Densmap,xslices,yslices,zslices,tblock,2*ftorder+1);
+ }
-if (bOut){
-outputfield(opt2fn("-o",NFILE,fnm),Densmap,xslices,yslices,zslices,tblock);
-}
+ if (bOut)
+ {
+ outputfield(opt2fn("-o",NFILE,fnm),Densmap,xslices,yslices,zslices,tblock);
+ }
-interfaces_txy(Densmap,xslices,yslices,zslices,tblock,binwz,eMeth,dens1,dens2,&surf1,&surf2,oenv);
+ interfaces_txy(Densmap,xslices,yslices,zslices,tblock,binwz,eMeth,dens1,dens2,&surf1,&surf2,oenv);
-if(bGraph){
+ if(bGraph)
+ {
-//Output surface-xpms
- nfxpm=opt2fns(&graphfiles,"-og",NFILE,fnm);
- if(nfxpm!=2){
- gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
- }
-writesurftoxpms(surf1, surf2, tblock,xslices,yslices,zslices,binw,binwz,graphfiles,zslices);
-}
+ /*Output surface-xpms*/
+ nfxpm=opt2fns(&graphfiles,"-og",NFILE,fnm);
+ if(nfxpm!=2)
+ {
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
+ }
+ writesurftoxpms(surf1, surf2, tblock,xslices,yslices,zslices,binw,binwz,graphfiles,zslices);
+ }
-//Output raw-data
-if (bRawOut){
- nfraw=opt2fns(&rawfiles,"-or",NFILE,fnm);
- if(nfraw!=2){
- gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
- }
- writeraw(surf1,surf2,tblock,xslices,yslices,rawfiles);
-}
+/*Output raw-data*/
+ if (bRawOut)
+ {
+ nfraw=opt2fns(&rawfiles,"-or",NFILE,fnm);
+ if(nfraw!=2)
+ {
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
+ }
+ writeraw(surf1,surf2,tblock,xslices,yslices,rawfiles);
+ }
-if(bFourier){
- nfspect=opt2fns(&spectra,"-Spect",NFILE,fnm);
- if(nfspect!=2){
- gmx_fatal(FARGS,"No or not correct number (2) of output-file-series: %d"
- ,nfspect);
- }
- powerspectavg_intf(surf1, surf2, tblock,xslices,yslices,spectra);
-}
+ if(bFourier)
+ {
+ nfspect=opt2fns(&spectra,"-Spect",NFILE,fnm);
+ if(nfspect!=2)
+ {
+ gmx_fatal(FARGS,"No or not correct number (2) of output-file-series: %d"
+ ,nfspect);
+ }
+ powerspectavg_intf(surf1, surf2, tblock,xslices,yslices,spectra);
+ }
-sfree(Densmap);
-if(bGraph || bFourier || bRawOut){
-sfree(surf1);
-sfree(surf2);
-}
+ sfree(Densmap);
+ if(bGraph || bFourier || bRawOut)
+ {
+ sfree(surf1);
+ sfree(surf2);
+ }
- thanx(stderr);
- return 0;
+ thanx(stderr);
+ return 0;
}
-/*
- * $Id: gmx_order.c,v 1.11 2009/03/11 12:13:15 spoel Exp $
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ * $Id: densorder.c,v 0.9
*
* This source code is part of
*
*
* GROningen MAchine for Chemical Simulations
*
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
+ * VERSION 3.0
+ *
+ * Copyright (c) 1991-2001
+ * BIOSON Research Institute, Dept. of Biophysical Chemistry
+ * University of Groningen, The Netherlands
+ *
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* To help us fund GROMACS development, we humbly ask that you cite
* the papers on the package - you can find them in the top README file.
*
- * For more info, check our website at http://www.gromacs.org
+ * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
*
* And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
+ * Gyas ROwers Mature At Cryogenic Speed
*/
+
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "matio.h"
#include "binsearch.h"
#include "powerspect.h"
-// #include "cmat.h"
-/****************************************************************************/
-/* This program calculates the order parameter per atom for an interface or */
-/* bilayer, averaged over time. */
-/* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
-/* It is assumed that the order parameter with respect to a box-axis */
-/* is calculated. In that case i = j = axis, and delta(ij) = 1. */
-/* */
-/* Peter Tieleman, April 1995 */
-/* P.J. van Maaren, November 2005 Added tetrahedral stuff */
-/****************************************************************************/
-
-
-
-
-
/* Print name of first atom in all groups in index file */
static void print_types(atom_id index[], atom_id a[], int ngrps,
- char *groups[], t_topology *top)
+ char *groups[], t_topology *top)
{
- int i;
+ int i;
- fprintf(stderr,"Using following groups: \n");
- for(i = 0; i < ngrps; i++)
- fprintf(stderr,"Groupname: %s First atomname: %s First atomnr %u\n",
- groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
- fprintf(stderr,"\n");
+ fprintf(stderr,"Using following groups: \n");
+ for(i = 0; i < ngrps; i++)
+ fprintf(stderr,"Groupname: %s First atomname: %s First atomnr %u\n",
+ groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
+ fprintf(stderr,"\n");
}
static void check_length(real length, int a, int b)
{
- if (length > 0.3)
- fprintf(stderr,"WARNING: distance between atoms %d and "
- "%d > 0.3 nm (%f). Index file might be corrupt.\n",
- a, b, length);
+ if (length > 0.3)
+ fprintf(stderr,"WARNING: distance between atoms %d and "
+ "%d > 0.3 nm (%f). Index file might be corrupt.\n",
+ a, b, length);
}
static void find_tetra_order_grid(t_topology top, int ePBC,
- int natoms, matrix box,
- rvec x[],int maxidx,atom_id index[],
- real time, real *sgmean, real *skmean,
- int nslicex, int nslicey, int nslicez,
- real ***sggrid,real ***skgrid)
+ int natoms, matrix box,
+ rvec x[],int maxidx,atom_id index[],
+ real time, real *sgmean, real *skmean,
+ int nslicex, int nslicey, int nslicez,
+ real ***sggrid,real ***skgrid)
{
- int ix,jx,i,j,k,l,n,*nn[4];
- rvec dx,rj,rk,urk,urj;
- real cost,cost2,*sgmol,*skmol,rmean,rmean2,r2,box2,*r_nn[4];
- t_pbc pbc;
- //t_mat *dmat;
- //t_dist *d;
- int slindex_x,slindex_y,slindex_z;
- int ***sl_count;
- real onethird=1.0/3.0;
- gmx_rmpbc_t gpbc;
+ int ix,jx,i,j,k,l,n,*nn[4];
+ rvec dx,rj,rk,urk,urj;
+ real cost,cost2,*sgmol,*skmol,rmean,rmean2,r2,box2,*r_nn[4];
+ t_pbc pbc;
+ int slindex_x,slindex_y,slindex_z;
+ int ***sl_count;
+ real onethird=1.0/3.0;
+ gmx_rmpbc_t gpbc;
- /* dmat = init_mat(maxidx, FALSE); */
-
- box2 = box[XX][XX] * box[XX][XX];
-
- // Initialize expanded sl_count array
- snew(sl_count,nslicex);
- for(i=0;i<nslicex;i++){
- snew(sl_count[i],nslicey);
- for(j=0;j<nslicey;j++){
- snew(sl_count[i][j],nslicez);
- }
- }
+ /* dmat = init_mat(maxidx, FALSE); */
+
+ box2 = box[XX][XX] * box[XX][XX];
+
+ /* Initialize expanded sl_count array */
+ snew(sl_count,nslicex);
+ for(i=0;i<nslicex;i++)
+ {
+ snew(sl_count[i],nslicey);
+ for(j=0;j<nslicey;j++)
+ {
+ snew(sl_count[i][j],nslicez);
+ }
+ }
- for (i=0; (i<4); i++) {
- snew(r_nn[i],natoms);
- snew(nn[i],natoms);
+ for (i=0; (i<4); i++)
+ {
+ snew(r_nn[i],natoms);
+ snew(nn[i],natoms);
- for (j=0; (j<natoms); j++) {
- r_nn[i][j] = box2;
+ for (j=0; (j<natoms); j++)
+ {
+ r_nn[i][j] = box2;
+ }
}
- }
- snew(sgmol,maxidx);
- snew(skmol,maxidx);
-
- /* Must init pbc every step because of pressure coupling */
- gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
- gmx_rmpbc(gpbc,natoms,box,x);
-
- *sgmean = 0.0;
- *skmean = 0.0;
- l=0;
- for (i=0; (i<maxidx); i++) { /* loop over index file */
- ix = index[i];
- for (j=0; (j<maxidx); j++) {
- if (i == j) continue;
-
- jx = index[j];
+ snew(sgmol,maxidx);
+ snew(skmol,maxidx);
+
+ /* Must init pbc every step because of pressure coupling */
+ gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
+ gmx_rmpbc(gpbc,natoms,box,x);
+
+ *sgmean = 0.0;
+ *skmean = 0.0;
+ l=0;
+ for (i=0; (i<maxidx); i++)
+ { /* loop over index file */
+ ix = index[i];
+ for (j=0; (j<maxidx); j++)
+ {
+
+ if (i == j) continue;
+
+ jx = index[j];
- pbc_dx(&pbc,x[ix],x[jx],dx);
- r2=iprod(dx,dx);
-
- /* set_mat_entry(dmat,i,j,r2); */
-
- /* determine the nearest neighbours */
- if (r2 < r_nn[0][i]) {
- r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
- r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
- r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
- r_nn[0][i] = r2; nn[0][i] = j;
- } else if (r2 < r_nn[1][i]) {
- r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
- r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
- r_nn[1][i] = r2; nn[1][i] = j;
- } else if (r2 < r_nn[2][i]) {
- r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
- r_nn[2][i] = r2; nn[2][i] = j;
- } else if (r2 < r_nn[3][i]) {
- r_nn[3][i] = r2; nn[3][i] = j;
- }
- }
-
-
- /* calculate mean distance between nearest neighbours */
- rmean = 0;
- for (j=0; (j<4); j++) {
- r_nn[j][i] = sqrt(r_nn[j][i]);
- rmean += r_nn[j][i];
- }
- rmean /= 4;
+ pbc_dx(&pbc,x[ix],x[jx],dx);
+ r2=iprod(dx,dx);
+
+ /* set_mat_entry(dmat,i,j,r2); */
+
+ /* determine the nearest neighbours */
+ if (r2 < r_nn[0][i])
+ {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
+ r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
+ r_nn[0][i] = r2; nn[0][i] = j;
+ } else if (r2 < r_nn[1][i])
+ {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
+ r_nn[1][i] = r2; nn[1][i] = j;
+ } else if (r2 < r_nn[2][i])
+ {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r2; nn[2][i] = j;
+ } else if (r2 < r_nn[3][i])
+ {
+ r_nn[3][i] = r2; nn[3][i] = j;
+ }
+ }
+
+
+ /* calculate mean distance between nearest neighbours */
+ rmean = 0;
+ for (j=0; (j<4); j++)
+ {
+ r_nn[j][i] = sqrt(r_nn[j][i]);
+ rmean += r_nn[j][i];
+ }
+ rmean /= 4;
- n = 0;
- sgmol[i] = 0.0;
- skmol[i] = 0.0;
-
- /* Chau1998a eqn 3 */
- /* angular part tetrahedrality order parameter per atom */
- for (j=0; (j<3); j++) {
- for (k=j+1; (k<4); k++) {
- pbc_dx(&pbc,x[ix],x[index[nn[k][i]]],rk);
- pbc_dx(&pbc,x[ix],x[index[nn[j][i]]],rj);
-
- unitv(rk,urk);
- unitv(rj,urj);
+ n = 0;
+ sgmol[i] = 0.0;
+ skmol[i] = 0.0;
+
+ /* Chau1998a eqn 3 */
+ /* angular part tetrahedrality order parameter per atom */
+ for (j=0; (j<3); j++)
+ {
+ for (k=j+1; (k<4); k++)
+ {
+ pbc_dx(&pbc,x[ix],x[index[nn[k][i]]],rk);
+ pbc_dx(&pbc,x[ix],x[index[nn[j][i]]],rj);
+
+ unitv(rk,urk);
+ unitv(rj,urj);
- cost = iprod(urk,urj) + onethird;
- cost2 = cost * cost;
-
- sgmol[i] += cost2;
- l++;
- n++;
- }
- }
- /* normalize sgmol between 0.0 and 1.0 */
- sgmol[i] = 3*sgmol[i]/32;
- *sgmean += sgmol[i];
-
- /* distance part tetrahedrality order parameter per atom */
- rmean2 = 4 * 3 * rmean * rmean;
- for (j=0; (j<4); j++) {
- skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
- /* printf("%d %f (%f %f %f %f) \n",
- i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
- */
- }
+ cost = iprod(urk,urj) + onethird;
+ cost2 = cost * cost;
+
+ sgmol[i] += cost2;
+ l++;
+ n++;
+ }
+ }
+ /* normalize sgmol between 0.0 and 1.0 */
+ sgmol[i] = 3*sgmol[i]/32;
+ *sgmean += sgmol[i];
+
+ /* distance part tetrahedrality order parameter per atom */
+ rmean2 = 4 * 3 * rmean * rmean;
+ for (j=0; (j<4); j++)
+ {
+ skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
+ /* printf("%d %f (%f %f %f %f) \n",
+ i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
+ */
+ }
- *skmean += skmol[i];
+ *skmean += skmol[i];
- /* Compute sliced stuff in x y z*/
- slindex_x = gmx_nint((1+x[i][XX]/box[XX][XX])*nslicex) % nslicex;
- slindex_y = gmx_nint((1+x[i][YY]/box[YY][YY])*nslicey) % nslicey;
- slindex_z = gmx_nint((1+x[i][ZZ]/box[ZZ][ZZ])*nslicez) % nslicez;
- sggrid[slindex_x][slindex_y][slindex_z] += sgmol[i];
- skgrid[slindex_x][slindex_y][slindex_z] += skmol[i];
- (sl_count[slindex_x][slindex_y][slindex_z])++;
- } /* loop over entries in index file */
+ /* Compute sliced stuff in x y z*/
+ slindex_x = gmx_nint((1+x[i][XX]/box[XX][XX])*nslicex) % nslicex;
+ slindex_y = gmx_nint((1+x[i][YY]/box[YY][YY])*nslicey) % nslicey;
+ slindex_z = gmx_nint((1+x[i][ZZ]/box[ZZ][ZZ])*nslicez) % nslicez;
+ sggrid[slindex_x][slindex_y][slindex_z] += sgmol[i];
+ skgrid[slindex_x][slindex_y][slindex_z] += skmol[i];
+ (sl_count[slindex_x][slindex_y][slindex_z])++;
+ } /* loop over entries in index file */
- *sgmean /= maxidx;
- *skmean /= maxidx;
+ *sgmean /= maxidx;
+ *skmean /= maxidx;
- for(i=0; (i<nslicex); i++) {
- for(j=0; j<nslicey; j++){
- for(k=0;k<nslicez;k++){
- if (sl_count[i][j][k] > 0) {
+ for(i=0; (i<nslicex); i++)
+ {
+ for(j=0; j<nslicey; j++)
+ {
+ for(k=0;k<nslicez;k++)
+ {
+ if (sl_count[i][j][k] > 0)
+ {
sggrid[i][j][k] /= sl_count[i][j][k];
skgrid[i][j][k] /= sl_count[i][j][k];
}
- }
- }
- }
-
- sfree(sl_count);
- sfree(sgmol);
- sfree(skmol);
- for (i=0; (i<4); i++) {
- sfree(r_nn[i]);
- sfree(nn[i]);
- }
+ }
+ }
+ }
+
+ sfree(sl_count);
+ sfree(sgmol);
+ sfree(skmol);
+ for (i=0; (i<4); i++)
+ {
+ sfree(r_nn[i]);
+ sfree(nn[i]);
+ }
}
/*Determines interface from tetrahedral order parameter in box with specified binwidth. */
/*Outputs interface positions(bins), the number of timeframes, and the number of surface-mesh points in xy*/
static void calc_tetra_order_interface(const char *fnNDX,const char *fnTPS,const char *fnTRX, real binw, int tblock,
- int *nframes, int *nslicex, int *nslicey,
- real sgang1,real sgang2,real ****intfpos,
- output_env_t oenv)
+ int *nframes, int *nslicex, int *nslicey,
+ real sgang1,real sgang2,real ****intfpos,
+ output_env_t oenv)
{
- FILE *fpsg=NULL,*fpsk=NULL;
- char *sgslfn="sg_ang_mesh"; //Hardcoded filenames for debugging
- char *skslfn="sk_dist_mesh";
- t_topology top;
- int ePBC;
- char title[STRLEN],subtitle[STRLEN];
- t_trxstatus *status;
- int natoms;
- real t;
- rvec *xtop,*x;
- matrix box;
- real sg,sk, sgintf, pos;
- atom_id **index;
- char **grpname;
- int i,j,k,n,*isize,ng, nslicez, framenr;
- real ***sg_grid,***sk_grid, ***sg_fravg, ***sk_fravg, ****sk_4d, ****sg_4d;
- int *perm;
- int ndx1, ndx2;
- int bins;
- const real onehalf=1.0/2.0;
- /* real ***intfpos[2]; pointers to arrays of two interface positions zcoord(framenr,xbin,ybin): intfpos[interface_index][t][nslicey*x+y]
- * i.e 1D Row-major order in (t,x,y) */
+ FILE *fpsg=NULL,*fpsk=NULL;
+ char *sgslfn="sg_ang_mesh"; /* Hardcoded filenames for debugging*/
+ char *skslfn="sk_dist_mesh";
+ t_topology top;
+ int ePBC;
+ char title[STRLEN],subtitle[STRLEN];
+ t_trxstatus *status;
+ int natoms;
+ real t;
+ rvec *xtop,*x;
+ matrix box;
+ real sg,sk, sgintf, pos;
+ atom_id **index;
+ char **grpname;
+ int i,j,k,n,*isize,ng, nslicez, framenr;
+ real ***sg_grid,***sk_grid, ***sg_fravg, ***sk_fravg, ****sk_4d, ****sg_4d;
+ int *perm;
+ int ndx1, ndx2;
+ int bins;
+ const real onehalf=1.0/2.0;
+ /* real ***intfpos[2]; pointers to arrays of two interface positions zcoord(framenr,xbin,ybin): intfpos[interface_index][t][nslicey*x+y]
+ * i.e 1D Row-major order in (t,x,y) */
- read_tps_conf(fnTPS,title,&top,&ePBC,&xtop,NULL,box,FALSE);
+ read_tps_conf(fnTPS,title,&top,&ePBC,&xtop,NULL,box,FALSE);
- *nslicex= (int)(box[XX][XX]/binw + onehalf); //Calculate slicenr from binwidth
- *nslicey= (int)(box[YY][YY]/binw + onehalf);
- nslicez= (int)(box[ZZ][ZZ]/binw + onehalf);
+ *nslicex= (int)(box[XX][XX]/binw + onehalf); /*Calculate slicenr from binwidth*/
+ *nslicey= (int)(box[YY][YY]/binw + onehalf);
+ nslicez= (int)(box[ZZ][ZZ]/binw + onehalf);
- ng = 1;
- /* get index groups */
- printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
- snew(grpname,ng);
- snew(index,ng);
- snew(isize,ng);
- get_index(&top.atoms,fnNDX,ng,isize,index,grpname);
-
- /* Analyze trajectory */
- natoms=read_first_x(oenv,&status,fnTRX,&t,&x,box);
- if ( natoms > top.atoms.nr )
- gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",
- top.atoms.nr,natoms);
- check_index(NULL,ng,index[0],NULL,natoms);
+ ng = 1;
+ /* get index groups */
+ printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
+ snew(grpname,ng);
+ snew(index,ng);
+ snew(isize,ng);
+ get_index(&top.atoms,fnNDX,ng,isize,index,grpname);
+
+ /* Analyze trajectory */
+ natoms=read_first_x(oenv,&status,fnTRX,&t,&x,box);
+ if ( natoms > top.atoms.nr )
+ gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",
+ top.atoms.nr,natoms);
+ check_index(NULL,ng,index[0],NULL,natoms);
/*Prepare structures for temporary storage of frame info*/
- snew(sg_grid,*nslicex);
- snew(sk_grid,*nslicex);
- for(i=0;i<*nslicex;i++){
- snew(sg_grid[i],*nslicey);
- snew(sk_grid[i],*nslicey);
- for(j=0;j<*nslicey;j++){
- snew(sg_grid[i][j],nslicez);
- snew(sk_grid[i][j],nslicez);
- }
- }
+ snew(sg_grid,*nslicex);
+ snew(sk_grid,*nslicex);
+ for(i=0;i<*nslicex;i++)
+ {
+ snew(sg_grid[i],*nslicey);
+ snew(sk_grid[i],*nslicey);
+ for(j=0;j<*nslicey;j++)
+ {
+ snew(sg_grid[i][j],nslicez);
+ snew(sk_grid[i][j],nslicez);
+ }
+ }
- sg_4d=NULL;
- sk_4d=NULL;
- *nframes = 0;
- framenr=0;
+ sg_4d=NULL;
+ sk_4d=NULL;
+ *nframes = 0;
+ framenr=0;
/* Loop over frames*/
- do {
-
-//Initialize box meshes (temporary storage for each tblock frame -reinitialise every tblock steps
- if(framenr%tblock ==0){
- srenew(sk_4d,*nframes+1);
- srenew(sg_4d,*nframes+1);
- snew(sg_fravg,*nslicex);
- snew(sk_fravg,*nslicex);
- for(i=0;i<*nslicex;i++){
- snew(sg_fravg[i],*nslicey);
- snew(sk_fravg[i],*nslicey);
- for(j=0;j<*nslicey;j++){
- snew(sg_fravg[i][j],nslicez);
- snew(sk_fravg[i][j],nslicez);
- }
- }
- }
+ do
+ {
+ /*Initialize box meshes (temporary storage for each tblock frame -reinitialise every tblock steps */
+ if(framenr%tblock ==0)
+ {
+ srenew(sk_4d,*nframes+1);
+ srenew(sg_4d,*nframes+1);
+ snew(sg_fravg,*nslicex);
+ snew(sk_fravg,*nslicex);
+ for(i=0;i<*nslicex;i++)
+ {
+ snew(sg_fravg[i],*nslicey);
+ snew(sk_fravg[i],*nslicey);
+ for(j=0;j<*nslicey;j++)
+ {
+ snew(sg_fravg[i][j],nslicez);
+ snew(sk_fravg[i][j],nslicez);
+ }
+ }
+ }
- find_tetra_order_grid(top,ePBC,natoms,box,x,isize[0],index[0],t,
- &sg,&sk,*nslicex,*nslicey,nslicez,sg_grid,sk_grid);
- for(i=0;i<*nslicex;i++){
- for(j=0;j<*nslicey;j++){
- for(k=0;k<nslicez;k++){
- sk_fravg[i][j][k]+=sk_grid[i][j][k]/tblock;
- sg_fravg[i][j][k]+=sg_grid[i][j][k]/tblock;
- }
- }
- }
-
- framenr++;
-
- if(framenr%tblock== 0){
- sk_4d[*nframes] = sk_fravg;
- sg_4d[*nframes] = sg_fravg;
+ find_tetra_order_grid(top,ePBC,natoms,box,x,isize[0],index[0],t,
+ &sg,&sk,*nslicex,*nslicey,nslicez,sg_grid,sk_grid);
+ for(i=0;i<*nslicex;i++)
+ {
+ for(j=0;j<*nslicey;j++)
+ {
+ for(k=0;k<nslicez;k++)
+ {
+ sk_fravg[i][j][k]+=sk_grid[i][j][k]/tblock;
+ sg_fravg[i][j][k]+=sg_grid[i][j][k]/tblock;
+ }
+ }
+ }
+
+ framenr++;
+
+ if(framenr%tblock== 0)
+ {
+ sk_4d[*nframes] = sk_fravg;
+ sg_4d[*nframes] = sg_fravg;
(*nframes)++;
}
- } while (read_next_x(oenv,status,&t,natoms,x,box));
-close_trj(status);
+ } while (read_next_x(oenv,status,&t,natoms,x,box));
+ close_trj(status);
- sfree(grpname);
- sfree(index);
- sfree(isize);
-
-//Debugging for printing out the entire order parameter meshes.
- if(debug){
- fpsg = xvgropen(sgslfn,"S\\sg\\N Angle Order Parameter / Meshpoint","(nm)","S\\sg\\N",oenv);
- fpsk = xvgropen(skslfn,"S\\sk\\N Distance Order Parameter / Meshpoint","(nm)","S\\sk\\N",oenv);
- for(n=0;n<(*nframes);n++){
- fprintf(fpsg,"%i\n",n);
- fprintf(fpsk,"%i\n",n);
- for(i=0; (i<*nslicex); i++) {
- for(j=0;j<*nslicey;j++){
- for(k=0;k<nslicez;k++){
- fprintf(fpsg,"%4f %4f %4f %8f\n",(i+0.5)*box[XX][XX]/(*nslicex),(j+0.5)*box[YY][YY]/(*nslicey),(k+0.5)*box[ZZ][ZZ]/nslicez,sg_4d[n][i][j][k]);
- fprintf(fpsk,"%4f %4f %4f %8f\n",(i+0.5)*box[XX][XX]/(*nslicex),(j+0.5)*box[YY][YY]/(*nslicey),(k+0.5)*box[ZZ][ZZ]/nslicez,sk_4d[n][i][j][k]);
- }
- }
- }
- }
-fclose(fpsg);
-fclose(fpsk);
- }
+ sfree(grpname);
+ sfree(index);
+ sfree(isize);
+
+ /*Debugging for printing out the entire order parameter meshes.*/
+ if(debug)
+ {
+ fpsg = xvgropen(sgslfn,"S\\sg\\N Angle Order Parameter / Meshpoint","(nm)","S\\sg\\N",oenv);
+ fpsk = xvgropen(skslfn,"S\\sk\\N Distance Order Parameter / Meshpoint","(nm)","S\\sk\\N",oenv);
+ for(n=0;n<(*nframes);n++)
+ {
+ fprintf(fpsg,"%i\n",n);
+ fprintf(fpsk,"%i\n",n);
+ for(i=0; (i<*nslicex); i++)
+ {
+ for(j=0;j<*nslicey;j++)
+ {
+ for(k=0;k<nslicez;k++)
+ {
+ fprintf(fpsg,"%4f %4f %4f %8f\n",(i+0.5)*box[XX][XX]/(*nslicex),(j+0.5)*box[YY][YY]/(*nslicey),(k+0.5)*box[ZZ][ZZ]/nslicez,sg_4d[n][i][j][k]);
+ fprintf(fpsk,"%4f %4f %4f %8f\n",(i+0.5)*box[XX][XX]/(*nslicex),(j+0.5)*box[YY][YY]/(*nslicey),(k+0.5)*box[ZZ][ZZ]/nslicez,sk_4d[n][i][j][k]);
+ }
+ }
+ }
+ }
+ fclose(fpsg);
+ fclose(fpsk);
+ }
- /* Find positions of interface z by scanning orderparam for each frame and for each xy-mesh cylinder along z*/
+ /* Find positions of interface z by scanning orderparam for each frame and for each xy-mesh cylinder along z*/
-//Simple trial: assume interface is in the middle of -sgang1 and sgang2
- sgintf=0.5*(sgang1+sgang2);
+ /*Simple trial: assume interface is in the middle of -sgang1 and sgang2*/
+ sgintf=0.5*(sgang1+sgang2);
-//Allocate memory for interface arrays;
-snew((*intfpos),2);
-snew((*intfpos)[0],*nframes);
-snew((*intfpos)[1],*nframes);
+ /*Allocate memory for interface arrays; */
+ snew((*intfpos),2);
+ snew((*intfpos)[0],*nframes);
+ snew((*intfpos)[1],*nframes);
-bins=(*nslicex)*(*nslicey);
+ bins=(*nslicex)*(*nslicey);
-snew(perm,nslicez); //permutation array for sorting along normal coordinate
+ snew(perm,nslicez); /*permutation array for sorting along normal coordinate*/
- for (n=0;n<*nframes;n++){
+ for (n=0;n<*nframes;n++)
+ {
snew((*intfpos)[0][n],bins);
snew((*intfpos)[1][n],bins);
- for(i=0;i<*nslicex;i++){
- for(j=0;j<*nslicey;j++){
- rangeArray(perm,nslicez); //reset permutation array to identity
- /*Binsearch returns 2 bin-numbers where the order param is <= setpoint sgintf*/
+ for(i=0;i<*nslicex;i++)
+ {
+ for(j=0;j<*nslicey;j++)
+ {
+ rangeArray(perm,nslicez); /*reset permutation array to identity*/
+ /*Binsearch returns 2 bin-numbers where the order param is <= setpoint sgintf*/
ndx1=start_binsearch(sg_4d[n][i][j],perm,0,nslicez/2-1,sgintf,1);
ndx2=start_binsearch(sg_4d[n][i][j],perm,nslicez/2,nslicez-1,sgintf,-1);
- /*Use linear interpolation to smooth out the interface position*/
+ /*Use linear interpolation to smooth out the interface position*/
- //left interface (0)
- /*if((sg_4d[n][i][j][perm[ndx1+1]]-sg_4d[n][i][j][perm[ndx1]])/sg_4d[n][i][j][perm[ndx1]] > 0.01){
- pos=( (sgintf-sg_4d[n][i][j][perm[ndx1]])*perm[ndx1+1]+(sg_4d[n][i][j][perm[ndx1+1]]-sgintf)*perm[ndx1 ])*/
- (*intfpos)[0][n][j+*nslicey*i]=(perm[ndx1]+onehalf)*binw;
- //right interface (1)
- //alpha=(sgintf-sg_4d[n][i][j][perm[ndx2]])/(sg_4d[n][i][j][perm[ndx2]+1]-sg_4d[n][i][j][perm[ndx2]]);
- //(*intfpos)[1][n][j+*nslicey*i]=((1-alpha)*perm[ndx2]+alpha*(perm[ndx2]+1)+onehalf)*box[ZZ][ZZ]/nslicez;
- (*intfpos)[1][n][j+*nslicey*i]=(perm[ndx2]+onehalf)*binw;
- }
- }
- }
+ /*left interface (0)*/
+ /*if((sg_4d[n][i][j][perm[ndx1+1]]-sg_4d[n][i][j][perm[ndx1]])/sg_4d[n][i][j][perm[ndx1]] > 0.01){
+ pos=( (sgintf-sg_4d[n][i][j][perm[ndx1]])*perm[ndx1+1]+(sg_4d[n][i][j][perm[ndx1+1]]-sgintf)*perm[ndx1 ])*/
+ (*intfpos)[0][n][j+*nslicey*i]=(perm[ndx1]+onehalf)*binw;
+ /*right interface (1)*/
+ /*alpha=(sgintf-sg_4d[n][i][j][perm[ndx2]])/(sg_4d[n][i][j][perm[ndx2]+1]-sg_4d[n][i][j][perm[ndx2]]);*/
+ /*(*intfpos)[1][n][j+*nslicey*i]=((1-alpha)*perm[ndx2]+alpha*(perm[ndx2]+1)+onehalf)*box[ZZ][ZZ]/nslicez;*/
+ (*intfpos)[1][n][j+*nslicey*i]=(perm[ndx2]+onehalf)*binw;
+ }
+ }
+ }
- //sfree(perm);
- sfree(sk_4d);
- sfree(sg_4d);
- //sfree(sg_grid);
- //sfree(sk_grid);
+ /*sfree(perm);*/
+ sfree(sk_4d);
+ sfree(sg_4d);
+ /*sfree(sg_grid);*/
+ /*sfree(sk_grid);*/
}
static void writesurftoxpms(real ***surf, int tblocks,int xbins, int ybins, real bw,char **outfiles,int maplevels )
{
- char numbuf[8];
- int n, i, j;
- real **profile1, **profile2;
- real max1, max2, min1, min2, *xticks, *yticks;
- t_rgb lo={1,1,1};
- t_rgb hi={0,0,0};
- FILE *xpmfile1, *xpmfile2;
+ char numbuf[8];
+ int n, i, j;
+ real **profile1, **profile2;
+ real max1, max2, min1, min2, *xticks, *yticks;
+ t_rgb lo={1,1,1};
+ t_rgb hi={0,0,0};
+ FILE *xpmfile1, *xpmfile2;
/*Prepare xpm structures for output*/
-//Allocate memory to tick's and matrices
-snew (xticks,xbins+1);
-snew (yticks,ybins+1);
+/*Allocate memory to tick's and matrices*/
+ snew (xticks,xbins+1);
+ snew (yticks,ybins+1);
-profile1=mk_matrix(xbins,ybins,FALSE);
-profile2=mk_matrix(xbins,ybins,FALSE);
+ profile1=mk_matrix(xbins,ybins,FALSE);
+ profile2=mk_matrix(xbins,ybins,FALSE);
-for (i=0;i<xbins+1;i++) xticks[i]+=bw;
-for (j=0;j<ybins+1;j++) yticks[j]+=bw;
-
-xpmfile1 = ffopen(outfiles[0],"w");
-xpmfile2 = ffopen(outfiles[1],"w");
-
-max1=max2=0.0;
-min1=min2=1000.00;
-
-for(n=0;n<tblocks;n++){
-sprintf(numbuf,"%5d",n);
-//Filling matrices for inclusion in xpm-files
- for(i=0;i<xbins;i++){
- for(j=0;j<ybins;j++){
+ for (i=0;i<xbins+1;i++) xticks[i]+=bw;
+ for (j=0;j<ybins+1;j++) yticks[j]+=bw;
+
+ xpmfile1 = ffopen(outfiles[0],"w");
+ xpmfile2 = ffopen(outfiles[1],"w");
+
+ max1=max2=0.0;
+ min1=min2=1000.00;
+
+ for(n=0;n<tblocks;n++)
+ {
+ sprintf(numbuf,"%5d",n);
+ /*Filling matrices for inclusion in xpm-files*/
+ for(i=0;i<xbins;i++)
+ {
+ for(j=0;j<ybins;j++)
+ {
profile1[i][j]=(surf[0][n][j+ybins*i]);
profile2[i][j]=(surf[1][n][j+ybins*i]);
- //Finding max and min values
+ /*Finding max and min values*/
if(profile1[i][j]>max1) max1=profile1[i][j];
if(profile1[i][j]<min1) min1=profile1[i][j];
if(profile2[i][j]>max2) max2=profile2[i][j];
if(profile2[i][j]<min2) min2=profile2[i][j];
- }
- }
+ }
+ }
- write_xpm(xpmfile1,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile1,min1,max1,lo,hi,&maplevels);
- write_xpm(xpmfile2,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile2,min2,max2,lo,hi,&maplevels);
-}
+ write_xpm(xpmfile1,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile1,min1,max1,lo,hi,&maplevels);
+ write_xpm(xpmfile2,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile2,min2,max2,lo,hi,&maplevels);
+ }
-ffclose(xpmfile1);
-ffclose(xpmfile2);
+ ffclose(xpmfile1);
+ ffclose(xpmfile2);
-sfree(profile1);
-sfree(profile2);
-sfree(xticks);
-sfree(yticks);
+ sfree(profile1);
+ sfree(profile2);
+ sfree(xticks);
+ sfree(yticks);
}
static void writeraw(real ***surf, int tblocks,int xbins, int ybins,char **fnms){
FILE *raw1, *raw2;
- int i,j,n;
+ int i,j,n;
raw1=ffopen(fnms[0],"w");
raw2=ffopen(fnms[1],"w");
fprintf(raw1,"#Legend\n#TBlock\n#Xbin Ybin Z t\n");
fprintf(raw2,"#Legend\n#TBlock\n#Xbin Ybin Z t\n");
- for (n=0;n<tblocks;n++){
+ for (n=0;n<tblocks;n++)
+ {
fprintf(raw1,"%5d\n",n);
fprintf(raw2,"%5d\n",n);
- for(i=0;i<xbins;i++){
- for(j=0;j<ybins;j++){
+ for(i=0;i<xbins;i++)
+ {
+ for(j=0;j<ybins;j++)
+ {
fprintf(raw1,"%i %i %8.5f\n",i,j,(surf[0][n][j+ybins*i]));
fprintf(raw2,"%i %i %8.5f\n",i,j,(surf[1][n][j+ybins*i]));
}
int gmx_hydorder(int argc,char *argv[])
{
- static const char *desc[] = {
- "The tetrahedrality order parameters can be determined",
- "around an atom. Both angle an distance order parameters are calculated. See",
- "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
- "for more details.[BR]"
- "This application calculates the orderparameter in a 3d-mesh in the box, and",
- "with 2 phases in the box gives the user the option to define a 2D interface in time",
- "separating the faces by specifying parameters -sgang1 and -sgang2 (It is important",
- "to select these judiciously)"};
+ static const char *desc[] = {
+ "The tetrahedrality order parameters can be determined",
+ "around an atom. Both angle an distance order parameters are calculated. See",
+ "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
+ "for more details.[BR]"
+ "This application calculates the orderparameter in a 3d-mesh in the box, and",
+ "with 2 phases in the box gives the user the option to define a 2D interface in time",
+ "separating the faces by specifying parameters -sgang1 and -sgang2 (It is important",
+ "to select these judiciously)"};
- int axis = 0;
- static int nsttblock=1;
- static int nlevels=100;
- static real binwidth = 1.0; /* binwidth in mesh */
- static real sg1 = 1;
- static real sg2 = 1; /* order parameters for bulk phases */
- static gmx_bool bFourier = FALSE;
- static gmx_bool bRawOut = FALSE;
- int frames,xslices,yslices; /* Dimensions of interface arrays*/
- real ***intfpos; /* Interface arrays (intfnr,t,xy) -potentially large */
- static char *normal_axis[] = { NULL, "z", "x", "y", NULL };
+ int axis = 0;
+ static int nsttblock=1;
+ static int nlevels=100;
+ static real binwidth = 1.0; /* binwidth in mesh */
+ static real sg1 = 1;
+ static real sg2 = 1; /* order parameters for bulk phases */
+ static gmx_bool bFourier = FALSE;
+ static gmx_bool bRawOut = FALSE;
+ int frames,xslices,yslices; /* Dimensions of interface arrays*/
+ real ***intfpos; /* Interface arrays (intfnr,t,xy) -potentially large */
+ static char *normal_axis[] = { NULL, "z", "x", "y", NULL };
- t_pargs pa[] = {
- { "-d", FALSE, etENUM, {normal_axis},
- "Direction of the normal on the membrane" },
- { "-bw", FALSE, etREAL, {&binwidth},
- "Binwidth of box mesh" },
- { "-sgang1",FALSE,etREAL, {&sg1},
- "tetrahedral angle parameter in Phase 1 (bulk)" },
- { "-sgang2",FALSE,etREAL, {&sg2},
- "tetrahedral angle parameter in Phase 2 (bulk)" },
- { "-tblock",FALSE,etINT,{&nsttblock},
- "Number of frames in one time-block average"},
- { "-nlevel", FALSE,etINT, {&nlevels},
- "Number of Height levels in 2D - XPixMaps"}
- };
-
- t_filenm fnm[] = { /* files for g_order */
- { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
- { efNDX, "-n", NULL, ffREAD }, /* index file */
- { efTPX, "-s", NULL, ffREAD }, /* topology file */
- { efXPM, "-o", "intf", ffWRMULT}, /* XPM- surface maps */
- { efOUT,"-or","raw", ffOPTWRMULT }, /* xvgr output file */
- { efOUT,"-Spect","intfspect",ffOPTWRMULT}, /* Fourier spectrum interfaces */
- };
+ t_pargs pa[] = {
+ { "-d", FALSE, etENUM, {normal_axis},
+ "Direction of the normal on the membrane" },
+ { "-bw", FALSE, etREAL, {&binwidth},
+ "Binwidth of box mesh" },
+ { "-sgang1",FALSE,etREAL, {&sg1},
+ "tetrahedral angle parameter in Phase 1 (bulk)" },
+ { "-sgang2",FALSE,etREAL, {&sg2},
+ "tetrahedral angle parameter in Phase 2 (bulk)" },
+ { "-tblock",FALSE,etINT,{&nsttblock},
+ "Number of frames in one time-block average"},
+ { "-nlevel", FALSE,etINT, {&nlevels},
+ "Number of Height levels in 2D - XPixMaps"}
+ };
+
+ t_filenm fnm[] = { /* files for g_order */
+ { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
+ { efNDX, "-n", NULL, ffREAD }, /* index file */
+ { efTPX, "-s", NULL, ffREAD }, /* topology file */
+ { efXPM, "-o", "intf", ffWRMULT}, /* XPM- surface maps */
+ { efOUT,"-or","raw", ffOPTWRMULT }, /* xvgr output file */
+ { efOUT,"-Spect","intfspect",ffOPTWRMULT}, /* Fourier spectrum interfaces */
+ };
#define NFILE asize(fnm)
- /*Filenames*/
- const char *ndxfnm,*tpsfnm,*trxfnm;
- char **spectra,**intfn, **raw;
- int nfspect,nfxpm, nfraw;
- output_env_t oenv;
+ /*Filenames*/
+ const char *ndxfnm,*tpsfnm,*trxfnm;
+ char **spectra,**intfn, **raw;
+ int nfspect,nfxpm, nfraw;
+ output_env_t oenv;
- CopyRight(stderr,argv[0]);
+ CopyRight(stderr,argv[0]);
- parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
- NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
- bFourier= opt2bSet("-Spect",NFILE,fnm);
- bRawOut = opt2bSet("-or",NFILE,fnm);
+ parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
+ NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
+ bFourier= opt2bSet("-Spect",NFILE,fnm);
+ bRawOut = opt2bSet("-or",NFILE,fnm);
- if (binwidth < 0.0)
- gmx_fatal(FARGS,"Can not have binwidth < 0");
+ if (binwidth < 0.0)
+ gmx_fatal(FARGS,"Can not have binwidth < 0");
- ndxfnm = ftp2fn(efNDX,NFILE,fnm);
- tpsfnm = ftp2fn(efTPX,NFILE,fnm);
- trxfnm = ftp2fn(efTRX,NFILE,fnm);
+ ndxfnm = ftp2fn(efNDX,NFILE,fnm);
+ tpsfnm = ftp2fn(efTPX,NFILE,fnm);
+ trxfnm = ftp2fn(efTRX,NFILE,fnm);
- /* Calculate axis */
- if (strcmp(normal_axis[0],"x") == 0) axis = XX;
- else if (strcmp(normal_axis[0],"y") == 0) axis = YY;
- else if (strcmp(normal_axis[0],"z") == 0) axis = ZZ;
- else gmx_fatal(FARGS,"Invalid axis, use x, y or z");
+ /* Calculate axis */
+ if (strcmp(normal_axis[0],"x") == 0) axis = XX;
+ else if (strcmp(normal_axis[0],"y") == 0) axis = YY;
+ else if (strcmp(normal_axis[0],"z") == 0) axis = ZZ;
+ else gmx_fatal(FARGS,"Invalid axis, use x, y or z");
- switch (axis) {
- case 0:
- fprintf(stderr,"Taking x axis as normal to the membrane\n");
- break;
- case 1:
- fprintf(stderr,"Taking y axis as normal to the membrane\n");
- break;
- case 2:
- fprintf(stderr,"Taking z axis as normal to the membrane\n");
- break;
- }
+ switch (axis) {
+ case 0:
+ fprintf(stderr,"Taking x axis as normal to the membrane\n");
+ break;
+ case 1:
+ fprintf(stderr,"Taking y axis as normal to the membrane\n");
+ break;
+ case 2:
+ fprintf(stderr,"Taking z axis as normal to the membrane\n");
+ break;
+ }
- /* tetraheder order parameter */
+ /* tetraheder order parameter */
/* If either of the options is set we compute both */
nfxpm=opt2fns(&intfn,"-o",NFILE,fnm);
- if(nfxpm!=2){
- gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
+ if(nfxpm!=2)
+ {
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
}
calc_tetra_order_interface(ndxfnm,tpsfnm,trxfnm,binwidth,nsttblock,&frames,&xslices,&yslices,sg1,sg2,&intfpos,oenv);
writesurftoxpms(intfpos,frames,xslices,yslices,binwidth,intfn,nlevels);
- if(bFourier){
- nfspect=opt2fns(&spectra,"-Spect",NFILE,fnm);
- if(nfspect!=2){
- gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfspect);
- }
- powerspectavg(intfpos, frames,xslices,yslices,spectra);
- }
+ if(bFourier)
+ {
+ nfspect=opt2fns(&spectra,"-Spect",NFILE,fnm);
+ if(nfspect!=2)
+ {
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfspect);
+ }
+ powerspectavg(intfpos, frames,xslices,yslices,spectra);
+ }
- if (bRawOut){
- nfraw=opt2fns(&raw,"-or",NFILE,fnm);
- if(nfraw!=2){
- gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfraw);
- }
- writeraw(intfpos,frames,xslices,yslices,raw);
- }
+ if (bRawOut)
+ {
+ nfraw=opt2fns(&raw,"-or",NFILE,fnm);
+ if(nfraw!=2)
+ {
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfraw);
+ }
+ writeraw(intfpos,frames,xslices,yslices,raw);
+ }
- thanx(stderr);
+ thanx(stderr);
- return 0;
+ return 0;
}
+/*
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Green Red Orange Magenta Azure Cyan Skyblue
+ */
+
#ifndef _interf_h
#define _interf_h
#include "typedefs.h"
typedef struct {
- real Z; /* Interface height-coordinate */
- real t; /*Interface thickness*/
+ real Z; /* Interface height-coordinate */
+ real t; /* Interface thickness */
} t_interf;
static void init_interf(t_interf *surf)
{
- surf->Z=0;
- surf->t=0;
+ surf->Z = 0;
+ surf->t = 0;
}
#endif
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ * $Id: densorder.c,v 0.9
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.0
+ *
+ * Copyright (c) 1991-2001
+ * BIOSON Research Institute, Dept. of Biophysical Chemistry
+ * University of Groningen, The Netherlands
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
+ *
+ * And Hey:
+ * Gyas ROwers Mature At Cryogenic Speed
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
/* Bi-directional linear search - iterative, returns the array dim of the closest element, assume array sorted small->large - also bidirectionally*/
#include <stdio.h>
-
+#include <types/simple.h>
int LinearSearch (double *array,int startindx, int stopindx, double key,int *count, int direction){
- //Iterative implementation - assume elements sorted
-int i;
-int keyindex;
-
-if(direction>=0){
- keyindex=startindx;
- for (i=startindx;i<=stopindx;i++){
- (*count)++;
- if(array[i]>key) {
- keyindex=i-1;
- return keyindex;
- }
+ /*Iterative implementation - assume elements sorted*/
+ int i;
+ int keyindex;
+
+ if(direction>=0){
+ keyindex=startindx;
+ for (i=startindx;i<=stopindx;i++){
+ (*count)++;
+ if(array[i]>key) {
+ keyindex=i-1;
+ return keyindex;
+ }
+ }
+ }
+ else if(direction<0 ){
+ keyindex=stopindx;
+ for(i=stopindx;i>=startindx;i--){
+ (*count)++;
+ if (array[i]>key){
+ keyindex=i+1;
+ return keyindex;
+ }
+ }
+ }
+
+ else printf("Error: startindex=stopindex");
+
+}
+
+
+/*Make range-array (Permutation identity) for sorting */
+void rangeArray(int *ar,int size)
+{
+ int i;
+ for (i=0;i<size;i++){
+ ar[i]=i;
}
}
-else if(direction<0 ){
- keyindex=stopindx;
- for(i=stopindx;i>=startindx;i--){
- (*count)++;
- if (array[i]>key){
- keyindex=i+1;
- return keyindex;
- }
+
+void pswap(int *v1, int *v2)
+{
+ int temp;
+ temp=*v1;
+ *v1=*v2;
+ *v2=temp;
+}
+
+
+void Swap (real *v1, real *v2)
+{
+ real temp;
+ temp = *v1;
+ *v1 = *v2;
+ *v2 = temp;
+}
+
+
+
+void insertionSort(real *arr, int *perm, int startndx, int endndx, int direction) {
+ int i, j;
+
+ if(direction>=0){
+ for (i = startndx; i <=endndx; i++) {
+ j = i;
+
+ while (j > startndx && arr[j - 1] > arr[j]) {
+ Swap(&arr[j],&arr[j-1]);
+ pswap(&perm[j],&perm[j-1]);
+ j--;
+ }
+
+ }
+
+ }
+
+ if(direction<0){
+ for (i = startndx; i <=endndx; i++) {
+ j = i;
+
+ while (j > startndx && arr[j - 1] < arr[j]) {
+ Swap(&arr[j],&arr[j-1]);
+ pswap(&perm[j],&perm[j-1]);
+ j--;
+ }
+
+ }
}
}
-else printf("Error: startindex=stopindex");
+int BinarySearch (real *array, int low, int high, real key,int direction)
+{
+ int mid, max, min;
+ max=high+2;
+ min=low+1;
+
+/*Iterative implementation*/
+
+ if (direction>=0){
+ while (max-min>1){
+ mid=(min+max)>>1;
+ if(key<array[mid-1]) max=mid;
+ else min=mid;
+ }
+ return min;
+ }
+
+ else if (direction<0){
+ while(max-min>1){
+ mid=(min+max)>>1;
+ if(key>array[mid-1]) max=mid;
+ else min=mid;
+ }
+ return min-1;
+
+ }/*end -ifelse direction*/
}
-
+
+int start_binsearch(real *array, int *perm, int low, int high,
+ real key, int direction)
+{
+ insertionSort(array,perm,low,high,direction);
+ return BinarySearch(array,low,high,key,direction);
+}
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ * $Id: densorder.c,v 0.9
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.0
+ *
+ * Copyright (c) 1991-2001
+ * BIOSON Research Institute, Dept. of Biophysical Chemistry
+ * University of Groningen, The Netherlands
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
+ *
+ * And Hey:
+ * Gyas ROwers Mature At Cryogenic Speed
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
#include <smalloc.h>
#include <gmx_fft.h>
#include <gmx_fatal.h>
void powerspectavg(real ***intftab, int tsteps, int xbins, int ybins, char **outfiles){
- //Fourier plans and output;
+ /*Fourier plans and output;*/
gmx_fft_t fftp;
t_complex *ftspect1; /* Spatial FFT of interface for each time frame and interface ftint[time,xycoord][0], ftintf[time,xycoord][1] for interface 1 and 2 respectively */
t_complex *ftspect2;
- real *pspectavg1; //power -spectrum 1st interface
- real *pspectavg2; // -------------- 2nd interface
+ real *pspectavg1; /*power -spectrum 1st interface*/
+ real *pspectavg2; /* -------------- 2nd interface*/
real *temp;
- FILE *datfile1,*datfile2; //data-files with interface data
- int n; //time index
- int fy=ybins/2+1; // number of (symmetric) fourier y elements;
- int rfl=xbins*fy; //length of real - DFT == Symmetric 2D matrix
+ FILE *datfile1,*datfile2; /*data-files with interface data*/
+ int n; /*time index*/
+ int fy=ybins/2+1; /* number of (symmetric) fourier y elements; */
+ int rfl=xbins*fy; /*length of real - DFT == Symmetric 2D matrix*/
int status;
/*Prepare data structures for FFT, with time averaging of power spectrum*/
gmx_fatal(status,__FILE__,__LINE__,"Error allocating FFT");
}
-//Initialize arrays
+/*Initialize arrays*/
snew(ftspect1,rfl);
snew(ftspect2,rfl);
snew(temp,rfl);
snew(pspectavg1,rfl);
snew(pspectavg2,rfl);
-//Fouriertransform directly (no normalization or anything)
-//NB! Check carefully indexes here
+/*Fouriertransform directly (no normalization or anything)*/
+/*NB! Check carefully indexes here*/
for (n=0;n<tsteps;n++){
gmx_fft_2d_real(fftp,GMX_FFT_REAL_TO_COMPLEX,intftab[0][n],ftspect1);
gmx_fft_2d_real(fftp,GMX_FFT_REAL_TO_COMPLEX,intftab[1][n],ftspect2);
- //Add to average for interface 1 here
+ /*Add to average for interface 1 here*/
addtoavgenergy(ftspect1,pspectavg1,rfl,tsteps);
- //Add to average for interface 2 here
+ /*Add to average for interface 2 here*/
addtoavgenergy(ftspect2,pspectavg2,rfl,tsteps);
}
- //Print out average energy-spectrum to outfiles[0] and outfiles[1];
+ /*Print out average energy-spectrum to outfiles[0] and outfiles[1];*/
datfile1 = ffopen(outfiles[0],"w");
datfile2 = ffopen(outfiles[1],"w");
-//Filling dat files with spectral data
+/*Filling dat files with spectral data*/
fprintf(datfile1,"%s\n","kx\t ky\t\tPower(kx,ky)");
fprintf(datfile2,"%s\n","kx\t ky\t\tPower(kx,ky)");
for(n=0;n<rfl;n++){
+/*
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Green Red Orange Magenta Azure Cyan Skyblue
+ */
+
#ifndef _powerspect_h
#define _powerspect_h
-#include <gmx_fft.h>
+#include "typedefs.h"
#include "interf.h"
-extern void powerspectavg(real ***interface, int t, int xbins, int ybins, char **outfiles);
-extern void powerspectavg_intf(t_interf ***if1,t_interf ***if2,int t, int xbins, int ybins, char **outfiles);
+extern void powerspectavg(real ***interface, int t, int xbins, int ybins,
+ char **outfiles);
+
+extern void powerspectavg_intf(t_interf ***if1,t_interf ***if2,int t,
+ int xbins, int ybins, char **outfiles);
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff