[gmx-users] g_rdf on Mac OS X 10.3, 10.4
David van der Spoel
spoel at xray.bmc.uu.se
Sat May 6 10:16:03 CEST 2006
Sukit Leekumjorn wrote:
> Hi
>
> I also have the exact same problem with g_rdf in Gromacs 3.3 and 3.3.1
> on OS X 10.3.9 (segmentation fault problem). Everything else seems to
> work fine. I was wandering if anyone encountered the same problem using
> the same OS and how to overcome it?
>
> Any help would be greatly appreciated
please submit a bugzilla with reproucible example
> Sukit
>
> Kay Gottschalk wrote:
>
>> I think we somehow fixed it in the attached version based on some
>> advise we found in the user list, but can't really remember the change
>> we made. Bad documentation. Try it out and see if it works.
>> Good luck,
>> Kay.
>>
>>
>> ------------------------------------------------------------------------
>>
>> /*
>> * $Id: gmx_rdf.c,v 1.1 2004/01/25 20:37:36 lindahl Exp $
>> * * 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
>> */
>> #ifdef HAVE_CONFIG_H
>> #include <config.h>
>> #endif
>>
>> #include <math.h>
>> #include <ctype.h>
>> #include "string2.h"
>> #include "sysstuff.h"
>> #include "typedefs.h"
>> #include "macros.h"
>> #include "vec.h"
>> #include "pbc.h"
>> #include "rmpbc.h"
>> #include "xvgr.h"
>> #include "copyrite.h"
>> #include "futil.h"
>> #include "statutil.h"
>> #include "tpxio.h"
>> #include "index.h"
>> #include "smalloc.h"
>> #include "fftgrid.h"
>> #include "calcgrid.h"
>> #include "nrnb.h"
>> #include "shift_util.h"
>> #include "pme.h"
>> #include "gstat.h"
>> #include "matio.h"
>>
>> typedef struct
>> {
>> const char *Label;
>> real a[4], b[4], c;
>> } t_CM_table;
>>
>> /*
>> * * f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
>> * i=1,4
>> */
>>
>> const t_CM_table CM_t[] =
>> {
>>
>> { "H", { 0.489918, 0.262003, 0.196767, 0.049879 },
>> { 20.6593, 7.74039, 49.5519, 2.20159 },
>> 0.001305 },
>> { "HO", { 0.489918, 0.262003, 0.196767, 0.049879 },
>> { 20.6593, 7.74039, 49.5519, 2.20159 },
>> 0.001305 },
>> { "HW", { 0.489918, 0.262003, 0.196767, 0.049879 },
>> { 20.6593, 7.74039, 49.5519, 2.20159 },
>> 0.001305 },
>> { "C", { 2.26069, 1.56165, 1.05075, 0.839259 },
>> { 22.6907, 0.656665, 9.75618, 55.5949 },
>> 0.286977 },
>> { "CB", { 2.26069, 1.56165, 1.05075, 0.839259 },
>> { 22.6907, 0.656665, 9.75618, 55.5949 },
>> 0.286977 },
>> { "CS1", { 0., 0., 0., 0.},{ 0., 0., 0., 0.}, 0.}, { "CS2", { 0.,
>> 0., 0., 0.},{ 0., 0., 0., 0.}, 0.},
>> { "N", { 12.2126, 3.13220, 2.01250, 1.16630 }, {
>> 0.005700, 9.89330, 28.9975, 0.582600 },
>> -11.529 },
>> { "O", { 3.04850, 2.28680, 1.54630, 0.867000 }, { 13.2771,
>> 5.70110, 0.323900, 32.9089 },
>> 0.250800 },
>> { "OW", { 3.04850, 2.28680, 1.54630, 0.867000 }, { 13.2771,
>> 5.70110, 0.323900, 32.9089 },
>> 0.250800 },
>> { "OWT3", { 3.04850, 2.28680, 1.54630, 0.867000 }, {
>> 13.2771, 5.70110, 0.323900, 32.9089 },
>> 0.250800 },
>> { "OA", { 3.04850, 2.28680, 1.54630, 0.867000 }, { 13.2771,
>> 5.70110, 0.323900, 32.9089 },
>> 0.250800 },
>> { "OS", { 3.04850, 2.28680, 1.54630, 0.867000 }, { 13.2771,
>> 5.70110, 0.323900, 32.9089 },
>> 0.250800 },
>> { "OSE", { 3.04850, 2.28680, 1.54630, 0.867000 },
>> { 13.2771, 5.70110, 0.323900, 32.9089 },
>> 0.250800 },
>> { "Na", { 3.25650, 3.93620, 1.39980, 1.00320 }, /* Na 1+ */
>> { 2.66710, 6.11530, 0.200100, 14.0390 }, 0.404000 },
>> { "CH3", { 0., 0., 0., 0.},{ 0., 0., 0., 0.}, 0.}, { "CH2", {
>> 0., 0., 0., 0.},{ 0., 0., 0., 0.}, 0.}, { "CH1", { 0., 0., 0.,
>> 0.},{ 0., 0., 0., 0.}, 0.}, };
>>
>> typedef struct
>> {
>> int n_angles;
>> int n_groups;
>> double lambda;
>> double energy;
>> double momentum;
>> double ref_k;
>> double **F;
>> int nSteps;
>> int total_n_atoms;
>> } structure_factor;
>>
>> typedef struct
>> {
>> rvec x;
>> int t;
>> } reduced_atom;
>>
>> real ** sf_table;
>>
>>
>> static void do_rdf(char *fnNDX,char *fnTPS,char *fnTRX,
>> char *fnRDF,char *fnCNRDF, char *fnHQ,
>> bool bCM,real cutoff,real binwidth,real fade)
>> {
>> FILE *fp;
>> int status;
>> char outf1[STRLEN],outf2[STRLEN];
>> char title[STRLEN];
>> int g,ng,natoms,i,j,k,nbin,j0,j1,n,nframes;
>> int **count;
>> char **grpname;
>> int *isize,isize_cm=0,nrdf=0,max_i;
>> atom_id **index,*index_cm=NULL;
>> unsigned long int *sum;
>> real t,boxmin,hbox,hbox2,cut2,r,r2,invbinw,normfac;
>> real segvol,spherevol,prev_spherevol,**rdf;
>> rvec *x,xcom,dx,*x_i1,xi;
>> real *inv_segvol,vol,vol_sum,rho;
>> bool *bExcl,bTop,bNonSelfExcl;
>> matrix box;
>> int **npairs;
>> atom_id ix,jx,***pairs;
>> t_topology top;
>> t_block *excl;
>> excl=NULL;
>> if (fnTPS) {
>> bTop=read_tps_conf(fnTPS,title,&top,&x,NULL,box,TRUE);
>> mk_single_top(&top);
>> if (bTop && !bCM)
>> /* get exclusions from topology */
>> excl=&(top.atoms.excl);
>> }
>> fprintf(stderr,"\nHow many groups do you want to calculate the RDF
>> of?\n");
>> do {
>> scanf("%d",&ng);
>> } while (ng < 1);
>> snew(grpname,ng+1);
>> snew(isize,ng+1);
>> snew(index,ng+1);
>> fprintf(stderr,"\nSelect a reference group and %d group%s\n",
>> ng,ng==1?"":"s");
>> if (fnTPS)
>> get_index(&top.atoms,fnNDX,ng+1,isize,index,grpname);
>> else
>> rd_index(fnNDX,ng+1,isize,index,grpname);
>> natoms=read_first_x(&status,fnTRX,&t,&x,box);
>> if ( !natoms )
>> fatal_error(0,"Could not read coordinates from statusfile\n");
>> if (fnTPS)
>> /* check with topology */
>> if ( natoms > top.atoms.nr ) fatal_error(0,"Trajectory (%d
>> atoms) does not match topology (%d atoms)",
>> natoms,top.atoms.nr);
>> /* check with index groups */
>> for (i=0; i<=ng; i++)
>> for (j=0; j<isize[i]; j++)
>> if ( index[i][j] >= natoms )
>> fatal_error(0,"Atom index (%d) in index group %s (%d atoms) larger "
>> "than number of atoms in trajectory (%d atoms)",
>> index[i][j],grpname[i],isize[i],natoms);
>> if (bCM) {
>> /* move index[0] to index_cm and make 'dummy' index[0] */
>> isize_cm=isize[0];
>> snew(index_cm,isize_cm);
>> for(i=0; i<isize[0]; i++)
>> index_cm[i]=index[0][i];
>> isize[0]=1;
>> index[0][0]=natoms;
>> srenew(index[0],isize[0]);
>> /* make space for center of mass */
>> srenew(x,natoms+1);
>> }
>> /* initialize some handy things */
>> boxmin = min( norm(box[XX]), min( norm(box[YY]), norm(box[ZZ]) ) );
>> hbox = boxmin / 2.0;
>> nbin = (int)(hbox / binwidth) + 1;
>> invbinw = 1.0 / binwidth;
>> hbox2 = sqr(hbox);
>> cut2 = sqr(cutoff);
>>
>> snew(count,ng);
>> snew(pairs,ng);
>> snew(npairs,ng);
>>
>> snew(bExcl,natoms);
>> max_i = 0;
>> for(g=0; g<ng; g++) {
>> if (isize[g+1] > max_i)
>> max_i = isize[g+1];
>>
>> /* this is THE array */
>> snew(count[g],nbin+1);
>> /* make pairlist array for groups and exclusions */
>> snew(pairs[g],isize[0]);
>> snew(npairs[g],isize[0]);
>> for(i=0; i<isize[0]; i++) {
>> ix = index[0][i];
>> for(j=0; j < natoms; j++)
>> bExcl[j] = FALSE;
>> /* exclusions? */
>> if (excl)
>> for( j = excl->index[ix]; j < excl->index[ix+1]; j++)
>> bExcl[excl->a[j]]=TRUE;
>> k = 0;
>> snew(pairs[g][i], isize[g+1]);
>> bNonSelfExcl = FALSE;
>> for(j=0; j<isize[g+1]; j++) {
>> jx = index[g+1][j];
>> if (!bExcl[jx])
>> pairs[g][i][k++]=jx;
>> else
>> /* Check if we have exclusions other than self exclusions */
>> bNonSelfExcl = bNonSelfExcl || (ix != jx);
>> }
>> if (bNonSelfExcl) {
>> npairs[g][i]=k;
>> srenew(pairs[g][i],npairs[g][i]);
>> } else {
>> /* Save a LOT of memory and some cpu cycles */
>> npairs[g][i]=-1;
>> sfree(pairs[g][i]);
>> }
>> }
>> }
>> sfree(bExcl);
>>
>> snew(x_i1,max_i);
>> nframes = 0;
>> vol_sum = 0;
>> do {
>> /* Must init pbc every step because of pressure coupling */
>> init_pbc(box);
>> rm_pbc(&top.idef,natoms,box,x,x);
>> vol = det(box);
>> vol_sum += vol;
>> if (bCM) {
>> /* calculate centre of mass */
>> clear_rvec(xcom);
>> for(i=0; (i < isize_cm); i++) {
>> ix = index_cm[i];
>> rvec_inc(xcom,x[ix]);
>> }
>> /* store it in the first 'group' */
>> for(j=0; (j<DIM); j++)
>> x[index[0][0]][j] = xcom[j] / isize_cm;
>> }
>>
>> for(g=0; g<ng; g++) {
>> /* Copy the indexed coordinates to a continuous array */
>> for(i=0; i<isize[g+1]; i++)
>> copy_rvec(x[index[g+1][i]],x_i1[i]);
>> for(i=0; i<isize[0]; i++) {
>> copy_rvec(x[index[0][i]],xi);
>> if (npairs[g][i] >= 0)
>> /* Expensive loop, because of indexing */
>> for(j=0; j<npairs[g][i]; j++) {
>> jx=pairs[g][i][j];
>> pbc_dx(xi,x[jx],dx);
>> r2=iprod(dx,dx);
>> if (r2>cut2 && r2<=hbox2)
>> count[g][(int)(sqrt(r2)*invbinw)]++;
>> }
>> else
>> /* Cheaper loop, no exclusions */
>> for(j=0; j<isize[g+1]; j++) {
>> pbc_dx(xi,x_i1[j],dx);
>> r2=iprod(dx,dx);
>> if (r2>cut2 && r2<=hbox2)
>> count[g][(int)(sqrt(r2)*invbinw)]++;
>> }
>> }
>> }
>> nframes++;
>> } while (read_next_x(status,&t,natoms,x,box));
>> fprintf(stderr,"\n");
>> close_trj(status);
>> sfree(x);
>> /* Average volume */
>> vol = vol_sum/nframes;
>> /* Calculate volume of sphere segments */
>> snew(inv_segvol,nbin);
>> prev_spherevol=0;
>> for(i=0; (i<nbin); i++) {
>> r = (i+1)*binwidth;
>> spherevol=(4.0/3.0)*M_PI*r*r*r;
>> segvol=spherevol-prev_spherevol;
>> inv_segvol[i]=1.0/segvol;
>> prev_spherevol=spherevol;
>> }
>> snew(rdf,ng);
>> for(g=0; g<ng; g++) {
>> /* We have to normalize by dividing by the number of frames */
>> rho = isize[g+1]/vol;
>> printf(""): /* DLB -- indicates a memory problem */
>> normfac = 1.0/((rho*nframes)*isize[0]);
>> /* Do the normalization */
>> nrdf = max(nbin-1,1+(2*fade/binwidth));
>> snew(rdf[g],nrdf);
>> for(i=0; (i<nbin-1); i++) {
>> r = (i+0.5)*binwidth;
>> if ((fade > 0) && (r >= fade))
>> rdf[g][i] =
>> 1+(count[g][i]*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
>> else
>> rdf[g][i] = count[g][i]*inv_segvol[i]*normfac;
>> }
>> for( ; (i<nrdf); i++)
>> rdf[g][i] = 1.0;
>> }
>> fp=xvgropen(fnRDF,"Radial Distribution","r","");
>> if (ng==1)
>> fprintf(fp,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
>> else {
>> fprintf(fp,"@ subtitle \"reference %s\"\n",grpname[0]);
>> xvgr_legend(fp,ng,grpname+1);
>> }
>> for(i=0; (i<nrdf); i++) {
>> fprintf(fp,"%10g", (i+0.5)*binwidth);
>> for(g=0; g<ng; g++)
>> fprintf(fp," %10g",rdf[g][i]);
>> fprintf(fp,"\n");
>> }
>> ffclose(fp);
>> do_view(fnRDF,NULL);
>>
>> /* h(Q) function: fourier transform of rdf */ if (fnHQ) {
>> int nhq = 401;
>> real *hq,*integrand,Q;
>> /* Get a better number density later! */
>> rho = isize[1]/vol;
>> snew(hq,nhq);
>> snew(integrand,nrdf);
>> for(i=0; (i<nhq); i++) {
>> Q = i*0.5;
>> integrand[0] = 0;
>> for(j=1; (j<nrdf); j++) {
>> r = (j+0.5)*binwidth;
>> integrand[j] = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
>> integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]-1.0);
>> }
>> hq[i] = print_and_integrate(debug,nrdf,binwidth,integrand,NULL,0);
>> }
>> fp=xvgropen(fnHQ,"h(Q)","Q(/nm)","h(Q)");
>> for(i=0; (i<nhq); i++) fprintf(fp,"%10g %10g\n",i*0.5,hq[i]);
>> ffclose(fp);
>> do_view(fnHQ,NULL);
>> sfree(hq);
>> sfree(integrand);
>> }
>> if (fnCNRDF) { normfac = 1.0/(isize[0]*nframes);
>> fp=xvgropen(fnCNRDF,"Cumulative Number RDF","r","number");
>> if (ng==1)
>> fprintf(fp,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
>> else {
>> fprintf(fp,"@ subtitle \"reference %s\"\n",grpname[0]);
>> xvgr_legend(fp,ng,grpname+1);
>> }
>> snew(sum,ng);
>> for(i=0; (i<nbin-1); i++) {
>> fprintf(fp,"%10g",i*binwidth);
>> for(g=0; g<ng; g++) {
>> fprintf(fp," %10g",(real)((double)sum[g]*normfac));
>> sum[g] += count[g][i];
>> }
>> fprintf(fp,"\n");
>> }
>> ffclose(fp);
>> sfree(sum);
>> do_view(fnCNRDF,NULL);
>> }
>>
>> for(g=0; g<ng; g++)
>> sfree(rdf[g]);
>> sfree(rdf);
>> }
>>
>> typedef struct {
>> int ndata;
>> real kkk;
>> real intensity;
>> } t_xdata;
>>
>> int comp_xdata(const void *a,const void *b)
>> {
>> t_xdata *xa,*xb;
>> real tmp;
>> xa = (t_xdata *)a;
>> xb = (t_xdata *)b;
>> if (xa->ndata == 0)
>> return 1;
>> else if (xb->ndata == 0)
>> return -1;
>> else {
>> tmp = xa->kkk - xb->kkk;
>> if (tmp < 0)
>> return -1;
>> else if (tmp > 0)
>> return 1;
>> else return 0;
>> }
>> }
>>
>> static t_xdata *init_xdata(int nx,int ny)
>> {
>> int ix,iy,i,j,maxkx,maxky;
>> t_xdata *data;
>> maxkx = (nx+1)/2;
>> maxky = (ny+1)/2;
>> snew(data,nx*ny);
>> for(i=0; (i<nx); i++) {
>> for(j=0; (j<ny); j++) {
>> ix = abs((i < maxkx) ? i : (i - nx)); iy = abs((j < maxky)
>> ? j : (j - ny)); data[ix*ny+iy].kkk = sqrt(ix*ix+iy*iy);
>> }
>> }
>> return data;
>> }
>>
>> static void extract_sq(t_fftgrid *fftgrid,int nbin,real k_max,real
>> lambda,
>> real count[],rvec box,int npixel,real *map[],
>> t_xdata data[])
>> {
>> int nx,ny,nz,nx2,ny2,nz2,la2,la12;
>> t_fft_c *ptr,*p0;
>> int i,j,k,maxkx,maxky,maxkz,n,ind,ix,iy;
>> real k1,kxy2,kz2,k2,z,kxy,kxy_max,cos_theta2,ttt,factor;
>> rvec lll,kk;
>> /*calc_lll(box,lll);
>> k_max = nbin/factor;
>> kxy_max = k_max/sqrt(3);*/
>> unpack_fftgrid(fftgrid,&nx,&ny,&nz,&nx2,&ny2,&nz2,
>> &la2,&la12,FALSE,(t_fft_r **)&ptr);
>> /* This bit copied from pme.c */
>> maxkx = (nx+1)/2;
>> maxky = (ny+1)/2;
>> maxkz = nz/2+1;
>> factor = nbin/k_max;
>> for(i=0; (i<nx); i++) {
>> #define IDX(i,n) (i<=n/2) ? (i) : (i-n)
>> kk[XX] = IDX(i,nx);
>> for(j=0; (j<ny); j++) {
>> kk[YY] = IDX(j,ny);
>> ind = INDEX(i,j,0);
>> p0 = ptr + ind;
>> for(k=0; (k<maxkz); k++,p0++) {
>> if ((i==0) && (j==0) && (k==0))
>> continue;
>> kk[ZZ] = IDX(k,nz);
>> z = sqrt(sqr(p0->re)+sqr(p0->im));
>> kxy2 = sqr(kk[XX]) + sqr(kk[YY]);
>> k2 = kxy2+sqr(kk[ZZ]);
>> k1 = sqrt(k2);
>> ind = k1*factor;
>> if (ind < nbin) {
>> /* According to:
>> * R. W. James (1962), * The Optical Principles of the
>> Diffraction of X-Rays,
>> * Oxbow press, Woodbridge Connecticut
>> * the intensity is proportional to (eq. 9.10):
>> * I = C (1+cos^2 [2 theta])/2 FFT
>> * And since
>> * cos[2 theta] = cos^2[theta] - sin^2[theta] = 2 cos^2[theta] -
>> 1 * we can compute the prefactor straight from cos[theta]
>> */
>> cos_theta2 = kxy2/k2;
>> /*ttt = z*0.5*(1+sqr(2*cos_theta2-1));*/
>> ttt = z*0.5*(1+cos_theta2);
>> count[ind] += ttt;
>> ix = ((i < maxkx) ? i : (i - nx)); iy = ((j < maxky) ? j :
>> (j - ny));
>> map[npixel/2+ix][npixel/2+iy] += ttt;
>> data[abs(ix)*ny+abs(iy)].ndata += 1;
>> data[abs(ix)*ny+abs(iy)].intensity += ttt;
>> }
>> else
>> fprintf(stderr,"k (%g) > k_max (%g)\n",k1,k_max);
>> }
>> }
>> }
>> }
>>
>> typedef struct {
>> char *name;
>> int nelec;
>> } t_element;
>>
>> static void do_sq(char *fnNDX,char *fnTPS,char *fnTRX,char *fnSQ,
>> char *fnXPM,real grid,real lambda,real distance,
>> int npixel,int nlevel)
>> {
>> FILE *fp;
>> t_element elem[] = { { "H", 1 }, { "C", 6 }, { "N", 7 }, { "O", 8
>> }, { "F", 9 }, { "S", 16 } };
>> #define NELEM asize(elem)
>> int status;
>> char title[STRLEN],*aname;
>> int natoms,i,j,k,nbin,j0,j1,n,nframes,pme_order;
>> real *count,**map;
>> char *grpname;
>> int isize;
>> atom_id *index;
>> real I0,C,t,k_max,factor,yfactor,segvol;
>> rvec *x,*xndx,box_size,kk,lll;
>> real fj0,*fj,max_spacing,r,lambda_1;
>> bool *bExcl;
>> matrix box;
>> int nx,ny,nz,nelectron;
>> atom_id ix,jx,**pairs;
>> splinevec *theta;
>> t_topology top;
>> t_fftgrid *fftgrid;
>> t_nrnb nrnb;
>> t_xdata *data;
>> /* bTop=read_tps_conf(fnTPS,title,&top,&x,NULL,box,TRUE); */
>>
>> fprintf(stderr,"\nSelect group for structure factor computation:\n");
>> get_index(&top.atoms,fnNDX,1,&isize,&index,&grpname);
>> natoms=read_first_x(&status,fnTRX,&t,&x,box);
>> if (isize < top.atoms.nr)
>> snew(xndx,isize);
>> else
>> xndx = x;
>> fprintf(stderr,"\n");
>> init_nrnb(&nrnb);
>> if ( !natoms )
>> fatal_error(0,"Could not read coordinates from statusfile\n");
>> /* check with topology */
>> if ( natoms > top.atoms.nr ) fatal_error(0,"Trajectory (%d
>> atoms) does not match topology (%d atoms)",
>> natoms,top.atoms.nr);
>>
>> /* Atomic scattering factors */
>> snew(fj,isize);
>> I0 = 0;
>> nelectron = 0;
>> for(i=0; (i<isize); i++) {
>> aname = *(top.atoms.atomname[index[i]]);
>> fj0 = 1;
>> if (top.atoms.atom[i].ptype == eptAtom) {
>> for(j=0; (j<NELEM); j++)
>> if (aname[0] == elem[j].name[0]) {
>> fj0 = elem[j].nelec;
>> break;
>> }
>> if (j == NELEM)
>> fprintf(stderr,"Warning: don't know number of electrons for atom
>> %s\n",aname);
>> }
>> /* Correct for partial charge */
>> fj[i] = fj0 - top.atoms.atom[index[i]].q;
>> nelectron += fj[i];
>> I0 += sqr(fj[i]);
>> }
>> if (debug) {
>> /* Dump scattering factors */
>> for(i=0; (i<isize); i++)
>> fprintf(debug,"Atom %3s-%5d q = %10.4f f = %10.4f\n",
>> *(top.atoms.atomname[index[i]]),index[i],
>> top.atoms.atom[index[i]].q,fj[i]);
>> }
>>
>> /* Constant for scattering */
>> C = sqr(1.0/(ELECTRONMASS_keV*KILO*ELECTRONVOLT*1e7*distance));
>> fprintf(stderr,"C is %g\n",C);
>> /* This bit is dimensionless */
>> nx = ny = nz = 0;
>> max_spacing = calc_grid(box,grid,&nx,&ny,&nz,1);
>> pme_order = max(4,1+(0.2/grid));
>> npixel = max(nx,ny);
>> data = init_xdata(nx,ny);
>> fprintf(stderr,"Largest grid spacing: %g nm, pme_order %d, %dx%d
>> pixel on image\n",
>> max_spacing,pme_order,npixel,npixel);
>> init_pme(stdout,NULL,nx,ny,nz,pme_order,isize,FALSE,eewg3D);
>> /* Determine largest k vector length. */
>> k_max = 1+sqrt(sqr(1+nx/2)+sqr(1+ny/2)+sqr(1+nz/2));
>>
>> /* this is the S(q) array */
>> nbin = npixel;
>> snew(count,nbin+1);
>> snew(map,npixel);
>> for(i=0; (i<npixel); i++)
>> snew(map[i],npixel);
>> nframes = 0;
>> do {
>> /* Put the atoms with scattering factor on a grid. Misuses
>> * an old routine from the PPPM code.
>> */
>> for(j=0; (j<DIM); j++)
>> box_size[j] = box[j][j];
>> /* Scale coordinates to the wavelength */
>> for(i=0; (i<isize); i++)
>> copy_rvec(x[index[i]],xndx[i]);
>> /* put local atoms on grid. */
>> fftgrid = spread_on_grid(stdout,isize,pme_order,xndx,fj,box,FALSE);
>>
>> /* FFT the density */
>> gmxfft3D(fftgrid,FFTW_FORWARD,NULL); /* Extract the Sq
>> function and sum it into the average array */
>> extract_sq(fftgrid,nbin,k_max,lambda,count,box_size,npixel,map,data);
>> nframes++;
>> } while (read_next_x(status,&t,natoms,x,box));
>> fprintf(stderr,"\n");
>> close_trj(status);
>> sfree(x);
>>
>> /* Normalize it ?? */ factor = k_max/(nbin);
>> yfactor = (1.0/nframes)/*(1.0/fftgrid->nxyz)*/;
>> fp=xvgropen(fnSQ,"Structure Factor","q (1/nm)","S(q)");
>> fprintf(fp,"@ subtitle \"Lambda = %g nm. Grid spacing = %g nm\"\n",
>> lambda,grid);
>> factor *= lambda;
>> for(i=0; i<nbin; i++) {
>> r = (i+0.5)*factor*2*M_PI;
>> segvol = 4*M_PI*sqr(r)*factor;
>> fprintf(fp,"%10g %10g\n",r,count[i]*yfactor/segvol);
>> }
>> ffclose(fp);
>> do_view(fnSQ,NULL);
>>
>> if (fnXPM) {
>> t_rgb rhi = { 0,0,0 }, rlo = { 1,1,1 };
>> real *tx,*ty,hi,inv_nframes;
>> hi = 0;
>> inv_nframes = 1.0/nframes;
>> snew(tx,npixel);
>> snew(ty,npixel);
>> for(i=0; (i<npixel); i++) {
>> tx[i] = i-npixel/2;
>> ty[i] = i-npixel/2;
>>
>> for(j=0; (j<npixel); j++) { map[i][j] *= inv_nframes;
>> hi = max(hi,map[i][j]);
>> }
>> }
>> fp = ffopen(fnXPM,"w");
>> write_xpm(fp,"Diffraction Image","Intensity","kx","ky",
>> nbin,nbin,tx,ty,map,0,hi,rlo,rhi,&nlevel);
>> fclose(fp);
>> sfree(tx);
>> sfree(ty);
>>
>> /* qsort(data,nx*ny,sizeof(data[0]),comp_xdata); fp =
>> ffopen("test.xvg","w");
>> for(i=0; (i<nx*ny); i++) {
>> if (data[i].ndata != 0) {
>> fprintf(fp,"%10.3f
>> %10.3f\n",data[i].kkk,data[i].intensity/data[i].ndata);
>> }
>> }
>> fclose(fp);
>> */
>> }
>> }
>>
>> t_complex *** rc_tensor_allocation(int x, int y, int z)
>> {
>> t_complex ***t;
>> int i,j;
>> snew(t,x);
>> t = (t_complex ***)calloc(x,sizeof(t_complex**));
>> if(!t) exit(fprintf(stderr,"\nallocation error"));
>> t[0] = (t_complex **)calloc(x*y,sizeof(t_complex*));
>> if(!t[0]) exit(fprintf(stderr,"\nallocation error"));
>> t[0][0] = (t_complex *)calloc(x*y*z,sizeof(t_complex));
>> if(!t[0][0]) exit(fprintf(stderr,"\nallocation error"));
>> for( j = 1 ; j < y ; j++) t[0][j] = t[0][j-1] + z;
>> for( i = 1 ; i < x ; i++) {
>> t[i] = t[i-1] + y;;
>> t[i][0] = t[i-1][0] + y*z;
>> for( j = 1 ; j < y ; j++) t[i][j] = t[i][j-1] + z;
>> }
>> return t;
>> }
>> int return_atom_type (char *type)
>> {
>> int i;
>> for (i = 0; (i < asize(CM_t)); i++)
>> if (!strcmp (type, CM_t[i].Label))
>> return i;
>> fatal_error(0,"\nError: atom type (%s) not in list (%d types
>> checked)!\n", type,i);
>>
>> return 0;
>> }
>>
>> double CMSF (int type, double lambda, double sin_theta)
>> /* * return Cromer-Mann fit for the atomic scattering factor:
>> * sin_theta is the sine of half the angle between incoming and scattered
>> * vectors. See g_sq.h for a short description of CM fit.
>> */
>> {
>> int i;
>> double tmp = 0.0, k2;
>>
>> k2 = (sqr (sin_theta) / sqr (10.0 * lambda));
>> /*
>> * united atoms case
>> * CH2 / CH3 groups */
>> if (!strcmp (CM_t[type].Label, "CS2") ||
>> !strcmp (CM_t[type].Label, "CH2") ||
>> !strcmp (CM_t[type].Label, "CH3") ||
>> !strcmp (CM_t[type].Label, "CS3")) {
>> tmp = CMSF (return_atom_type ("C"), lambda, sin_theta);
>> if (!strcmp (CM_t[type].Label, "CH3") ||
>> !strcmp (CM_t[type].Label, "CS3"))
>> tmp += 3.0 * CMSF (return_atom_type ("H"), lambda, sin_theta);
>> else
>> tmp += 2.0 * CMSF (return_atom_type ("H"), lambda, sin_theta);
>> }
>> /* all atom case */
>> else {
>> tmp = CM_t[type].c;
>> for (i = 0; i < 4; i++)
>> tmp += CM_t[type].a[i] * exp (-CM_t[type].b[i] * k2);
>> }
>> return tmp;
>> }
>>
>> void compute_scattering_factor_table (structure_factor * sf)
>> {
>> /*
>> * this function build up a table of scattering factors for every atom
>> * type and for every scattering angle.
>> */
>> int i, j;
>> double sin_theta,q;
>>
>> /* \hbar \omega \lambda = hc = 1239.842 eV * nm */
>> sf->momentum = ((double) (2. * 1000.0 * M_PI * sf->energy) /
>> 1239.842);
>> sf->lambda = 1239.842 / (1000.0 * sf->energy);
>> fprintf (stderr, "\nwavelenght = %f nm\n", sf->lambda);
>> snew (sf_table,asize (CM_t));
>> for (i = 0; (i < asize(CM_t)); i++) {
>> snew (sf_table[i], sf->n_angles);
>> for (j = 0; j < sf->n_angles; j++) {
>> q = ((double) j * sf->ref_k);
>> /* theta is half the angle between incoming and scattered wavevectors. */
>> sin_theta = q / (2.0 * sf->momentum);
>> sf_table[i][j] = CMSF (i, sf->lambda, sin_theta);
>> }
>> }
>> }
>>
>> int * create_indexed_atom_type (reduced_atom * atm, int size)
>> {
>> /* * create an index of the atom types found in a group
>> * i.e.: for water index_atp[0]=type_number_of_O and
>> * index_atp[1]=type_number_of_H
>> * * the last element is set to 0 */
>> int *index_atp, i, i_tmp, j;
>>
>> snew (index_atp, 1);
>> i_tmp = 1;
>> index_atp[0] = atm[0].t;
>> for (i = 1; i < size; i++) {
>> for (j = 0; j < i_tmp; j++)
>> if (atm[i].t == index_atp[j])
>> break;
>> if (j == i_tmp) { /* i.e. no indexed atom type is == to
>> atm[i].t */
>> i_tmp++;
>> srenew (index_atp, i_tmp * sizeof (int));
>> index_atp[i_tmp - 1] = atm[i].t;
>> }
>> }
>> i_tmp++;
>> srenew (index_atp, i_tmp * sizeof (int));
>> index_atp[i_tmp - 1] = 0;
>> return index_atp;
>> }
>>
>> void rearrange_atoms (reduced_atom * positions, t_trxframe * fr,
>> atom_id * index,
>> int isize, t_topology * top, real ** sf_table, bool flag)
>> /* given the group's index, return the (continuous) array of atoms */
>> {
>> int i;
>>
>> if (flag)
>> for (i = 0; i < isize; i++)
>> positions[i].t =
>> return_atom_type (*(top->atoms.atomtype[index[i]]));
>> for (i = 0; i < isize; i++)
>> copy_rvec (fr->x[index[i]], positions[i].x);
>> }
>>
>>
>> int atp_size (int *index_atp)
>> {
>> int i = 0;
>>
>> while (index_atp[i])
>> i++;
>> return i;
>> }
>>
>> void compute_structure_factor (structure_factor * sf, matrix box,
>> reduced_atom * red, int isize, real start_q,
>> real end_q, int group)
>> {
>> t_complex ***tmpSF;
>> rvec k_factor;
>> real kdotx, asf, kx, ky, kz, krr;
>> int kr, maxkx, maxky, maxkz, i, j, k, p, *counter;
>>
>>
>> k_factor[XX] = 2 * M_PI / box[XX][XX];
>> k_factor[YY] = 2 * M_PI / box[YY][YY];
>> k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ];
>>
>> maxkx = (int) rint (end_q / k_factor[XX]);
>> maxky = (int) rint (end_q / k_factor[YY]);
>> maxkz = (int) rint (end_q / k_factor[ZZ]);
>>
>> snew (counter, sf->n_angles);
>>
>> tmpSF = rc_tensor_allocation(maxkx,maxky,maxkz);
>> /*
>> * The big loop...
>> * compute real and imaginary part of the structure factor for every
>> * (kx,ky,kz)) */
>> fprintf(stderr,"\n");
>> for (i = 0; i < maxkx; i++) {
>> fprintf (stderr,"\rdone %3.1f%% ", (double)(100.0*(i+1))/maxkx);
>> fflush (stdout);
>> kx = i * k_factor[XX];
>> for (j = 0; j < maxky; j++) {
>> ky = j * k_factor[YY];
>> for (k = 0; k < maxkz; k++)
>> if (i != 0 || j != 0 || k != 0) {
>> kz = k * k_factor[ZZ];
>> krr = sqrt (sqr (kx) + sqr (ky) + sqr (kz));
>> if (krr >= start_q && krr <= end_q) {
>> kr = (int) rint (krr/sf->ref_k);
>> if (kr < sf->n_angles) {
>> counter[kr]++; /* will be used for the copmutation
>> of the average*/
>> for (p = 0; p < isize; p++) {
>> asf = sf_table[red[p].t][kr];
>>
>> kdotx = kx * red[p].x[XX] +
>> ky * red[p].x[YY] + kz * red[p].x[ZZ];
>>
>> tmpSF[i][j][k].re += cos (kdotx) * asf;
>> tmpSF[i][j][k].im += sin (kdotx) * asf;
>> }
>> }
>> }
>> }
>> }
>> } /* end loop on i */
>> /*
>> * compute the square modulus of the structure factor, averaging on
>> the surface
>> * kx*kx + ky*ky + kz*kz = krr*krr * note that this is correct only
>> for a (on the macroscopic scale)
>> * isotropic system. */
>> for (i = 0; i < maxkx; i++) {
>> kx = i * k_factor[XX]; for (j = 0; j < maxky; j++) {
>> ky = j * k_factor[YY]; for (k = 0; k < maxkz; k++) {
>> kz = k * k_factor[ZZ]; krr = sqrt (sqr (kx) + sqr (ky)
>> + sqr (kz)); if (krr >= start_q && krr <= end_q) {
>> kr = (int) rint (krr / sf->ref_k); if (kr <
>> sf->n_angles && counter[kr] != 0)
>> sf->F[group][kr] +=
>> (sqr (tmpSF[i][j][k].re) +
>> sqr (tmpSF[i][j][k].im))/ counter[kr];
>> }
>> }
>> }
>> } sfree (counter); free(tmpSF[0][0]); free(tmpSF[0]); free(tmpSF);
>> }
>>
>> void save_data (structure_factor * sf, char *file, int ngrps, real
>> start_q,
>> real end_q)
>> {
>>
>> FILE *fp;
>> int i, g = 0;
>> double *tmp, polarization_factor, A;
>>
>> fp = xvgropen (file, "Scattering Intensity", "q (1/nm)",
>> "Intensity (a.u.)");
>>
>> snew (tmp, ngrps);
>>
>> for (g = 0; g < ngrps; g++)
>> for (i = 0; i < sf->n_angles; i++) {
>> /*
>> * theta is half the angle between incoming and scattered
>> vectors.
>> * * polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 -
>> 0.5 * sin^2(2*theta)
>> * * sin(theta) = q/(2k) := A -> sin^2(theta) =
>> 4*A^2 (1-A^2) ->
>> * -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
>> */
>> A = (double) (i * sf->ref_k) / (2.0 * sf->momentum);
>> polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A));
>> sf->F[g][i] *= polarization_factor;
>> }
>> for (i = 0; i < sf->n_angles; i++) {
>> if (i * sf->ref_k >= start_q && i * sf->ref_k <= end_q) {
>> fprintf (fp, "%10.5f ", i * sf->ref_k);
>> for (g = 0; g < ngrps; g++)
>> fprintf (fp, " %10.5f ", (sf->F[g][i]) /(
>> sf->total_n_atoms*
>> sf->nSteps));
>> fprintf (fp, "\n");
>> }
>> }
>> ffclose (fp);
>> }
>>
>> int
>> do_scattering_intensity (char* fnTPS, char* fnNDX, char* fnXVG, char
>> *fnTRX,
>> real start_q,real end_q, real energy)
>> {
>> int i,ng,*isize,status,flags = TRX_READ_X,**index_atp;
>> char **grpname,title[STRLEN];
>> atom_id **index;
>> t_topology top;
>> t_trxframe fr;
>> reduced_atom **red;
>> structure_factor *sf;
>> rvec *xtop;
>> matrix box;
>> double r_tmp;
>>
>> snew (sf, 1);
>> sf->energy = energy;
>> /* Read the topology informations */
>> read_tps_conf (fnTPS, title, &top, &xtop, NULL, box, TRUE);
>> sfree (xtop);
>>
>> /* groups stuff... */
>> fprintf (stderr,
>> "\nHow many groups do you want to calculate the I(q) of?\n");
>> do {
>> scanf ("%d", &ng);
>> }
>> while (ng < 1);
>> snew (isize, ng);
>> snew (index, ng);
>> snew (grpname, ng);
>>
>> fprintf (stderr, "\nSelect %d group%s\n", ng,
>> ng == 1 ? "" : "s");
>> if (fnTPS)
>> get_index (&top.atoms, fnNDX, ng, isize, index, grpname);
>> else
>> rd_index (fnNDX, ng, isize, index, grpname);
>>
>> /* The first time we read data is a little special */
>> read_first_frame (&status, fnTRX, &fr, flags);
>>
>> sf->total_n_atoms = fr.natoms;
>> snew (red, ng);
>> snew (index_atp, ng);
>>
>> r_tmp = max (box[XX][XX], box[YY][YY]);
>> r_tmp = (double) max (box[ZZ][ZZ], r_tmp);
>>
>> sf->ref_k = (2.0 * M_PI) / (r_tmp);
>> /* ref_k will be the reference momentum unit */
>> sf->n_angles = (int) rint (end_q / sf->ref_k);
>>
>> snew (sf->F, ng);
>> for (i = 0; i < ng; i++)
>> snew (sf->F[i], sf->n_angles);
>> for (i = 0; i < ng; i++) {
>> snew (red[i], isize[i]);
>> rearrange_atoms (red[i], &fr, index[i], isize[i], &top, sf_table, 1);
>> index_atp[i] = create_indexed_atom_type (red[i], isize[i]);
>> }
>> compute_scattering_factor_table (sf);
>> /* This is the main loop over frames */
>>
>> do {
>> init_pbc (box);
>> sf->nSteps++;
>> for (i = 0; i < ng; i++) {
>> rearrange_atoms (red[i], &fr, index[i], isize[i], &top,
>> sf_table, 0);
>>
>> compute_structure_factor (sf, box, red[i], isize[i],
>> start_q, end_q, i);
>> }
>> }
>> while (read_next_frame (status, &fr));
>>
>> save_data (sf, fnXVG, ng, start_q, end_q);
>>
>> return 0;
>> }
>>
>> int gmx_rdf(int argc,char *argv[])
>> {
>> static char *desc[] = {
>> "The structure of liquids can be studied by either neutron or X-ray",
>> "scattering. The most common way to describe liquid structure is
>> by a",
>> "radial distribution function. However, this is not easy to obtain
>> from",
>> "a scattering experiment.[PAR]",
>> "g_rdf calculates radial distribution functions in different ways.",
>> "The normal method is around a (set of) particle(s), the other
>> method",
>> "is around the center of mass of a set of particles.[PAR]",
>> "If a run input file is supplied ([TT]-s[tt]), exclusions defined",
>> "in that file are taken into account when calculating the rdf.",
>> "The option [TT]-cut[tt] is meant as an alternative way to avoid",
>> "intramolecular peaks in the rdf plot.",
>> "It is however better to supply a run input file with a higher
>> number of",
>> "exclusions. For eg. benzene a topology with nrexcl set to 5",
>> "would eliminate all intramolecular contributions to the rdf.",
>> "Note that all atoms in the selected groups are used, also the ones",
>> "that don't have Lennard-Jones interactions.[PAR]",
>> "Option [TT]-cn[tt] produces the cumulative number rdf.[PAR]"
>> "To bridge the gap between theory and experiment structure factors
>> can",
>> "be computed (option [TT]-sq[tt]). The algorithm uses FFT, the grid"
>> "spacing of which is determined by option [TT]-grid[tt]."
>> };
>> static bool bCM=FALSE;
>> static real
>> cutoff=0,binwidth=0.001,grid=0.05,fade=0.0,lambda=0.1,distance=10;
>> static int npixel=256,nlevel=20;
>> static real start_q=0.0, end_q=60.0, energy=12.0;
>> t_pargs pa[] = {
>> { "-bin", FALSE, etREAL, {&binwidth},
>> "Binwidth (nm)" },
>> { "-com", FALSE, etBOOL, {&bCM},
>> "RDF with respect to the center of mass of first group" },
>> { "-cut", FALSE, etREAL, {&cutoff},
>> "Shortest distance (nm) to be considered"},
>> { "-fade", FALSE, etREAL, {&fade},
>> "From this distance onwards the RDF is tranformed by g'(r) = 1 +
>> [g(r)-1] exp(-(r/fade-1)^2 to make it go to 1 smoothly. If fade is 0.0
>> nothing is done." },
>> { "-grid", FALSE, etREAL, {&grid},
>> "[HIDDEN]Grid spacing (in nm) for FFTs when computing structure
>> factors" },
>> { "-npixel", FALSE, etINT, {&npixel},
>> "[HIDDEN]# pixels per edge of the square detector plate" },
>> { "-nlevel", FALSE, etINT, {&nlevel},
>> "Number of different colors in the diffraction image" },
>> { "-distance", FALSE, etREAL, {&distance},
>> "[HIDDEN]Distance (in cm) from the sample to the detector" },
>> { "-wave", FALSE, etREAL, {&lambda},
>> "[HIDDEN]Wavelength for X-rays/Neutrons for scattering. 0.1 nm
>> corresponds to roughly 12 keV" },
>> {"-startq", FALSE, etREAL, {&start_q},
>> "Starting q (1/nm) "},
>> {"-endq", FALSE, etREAL, {&end_q},
>> "Ending q (1/nm)"},
>> {"-energy", FALSE, etREAL, {&energy},
>> "Energy of the incoming X-ray (keV) "}
>> };
>> #define NPA asize(pa)
>> char *fnTPS,*fnNDX;
>> bool bSQ,bRDF;
>> t_filenm fnm[] = {
>> { efTRX, "-f", NULL, ffREAD },
>> { efTPS, NULL, NULL, ffOPTRD },
>> { efNDX, NULL, NULL, ffOPTRD },
>> { efXVG, "-o", "rdf", ffOPTWR },
>> { efXVG, "-sq", "sq", ffOPTWR },
>> { efXVG, "-cn", "rdf_cn", ffOPTWR },
>> { efXVG, "-hq", "hq", ffOPTWR },
>> /* { efXPM, "-image", "sq", ffOPTWR }*/
>> };
>> #define NFILE asize(fnm)
>> CopyRight(stderr,argv[0]);
>> parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
>> NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
>>
>> fnTPS = ftp2fn_null(efTPS,NFILE,fnm);
>> fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
>> /* bSQ = opt2bSet("-sq",NFILE,fnm) ||
>> opt2parg_bSet("-grid",NPA,pa); */
>> bSQ = opt2bSet("-sq",NFILE,fnm);
>> bRDF = opt2bSet("-o",NFILE,fnm) || !bSQ;
>> if (bSQ) {
>> if (!fnTPS)
>> fatal_error(0,"Need a tps file for calculating structure
>> factors\n");
>> }
>> else {
>> if (!fnTPS && !fnNDX)
>> fatal_error(0,"Neither index file nor topology file specified\n"
>> " Nothing to do!");
>> }
>>
>> if (bSQ)
>> do_scattering_intensity(fnTPS,fnNDX,opt2fn("-sq",NFILE,fnm),ftp2fn(efTRX,NFILE,fnm),
>>
>> start_q, end_q, energy );
>> /* old structure factor code */
>> /* do_sq(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),opt2fn("-sq",NFILE,fnm),
>> ftp2fn(efXPM,NFILE,fnm),grid,lambda,distance,npixel,nlevel);
>> */
>> if (bRDF) do_rdf(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),
>> opt2fn("-o",NFILE,fnm),opt2fn_null("-cn",NFILE,fnm),
>> opt2fn_null("-hq",NFILE,fnm),
>> bCM,cutoff,binwidth,fade);
>>
>> thanx(stderr);
>> return 0;
>> }
>>
>> ------------------------------------------------------------------------
>>
>>
>>
>>
>>
>> On May 4, 2006, at 9:50 PM, David van der Spoel wrote:
>>
>>> Jennifer Rendell wrote:
>>>
>>>> Dear friends,
>>>> I have seen a couple of postings (January 2004, gmx-developers,
>>>> David Bostick, and November 2003, gmx-users, Kay Gottschalk) on
>>>> segmentation faults from g_rdf on Mac OS X systems.
>>>> I am using the "Getting Started" section of the gromacs web pages,
>>>> in particular, the section on water, where the command g_rdf is called.
>>>> Below I show the results from two versions of gromacs, 3.2.1 on a
>>>> Mac OS X 10.4, and 3.1.5_pre2 on a Mac OS X 10.3. The results are
>>>> similar in that each results in a segmentation fault. g_rms and
>>>> gmxcheck work as I expect.
>>>> Any ideas on how to fix this? Jennifer
>>>
>>>
>>> Maybe you don't want to hear this, but please upgrade to 3.3.1 and
>>> try again. If the problem persists in the latest version please post
>>> a bugzilla.
>>>
>>>
>>>> *******************************************************************************
>>>> 1. gromacs 3.2.1 on Mac OS X 10.4 (Tiger). g_rdf gives the following
>>>> results:
>>>> % g_rdf -f water.trr -n oxygen.ndx -o rdf.xvg -s water.tpr
>>>> :-) G R O M A C S (-:
>>>> Giving Russians Opium May Alter Current Situation
>>>> :-) VERSION 3.2.1 (-:
>>>> (deleted some lines)
>>>> :-) g_rdf (-:
>>>> Option Filename Type Description
>>>> ------------------------------------------------------------
>>>> -f water.trr Input Generic trajectory: xtc trr trj
>>>> gro g96 pdb
>>>> -s water.tpr Input, Opt! Structure+mass(db): tpr tpb tpa
>>>> gro g96 pdb
>>>> xml
>>>> -n oxygen.ndx Input, Opt! Index file
>>>> -o rdf.xvg Output, Opt! xvgr/xmgr file
>>>> (deleted all remaining option lines)
>>>> Reading file water.tpr, VERSION 3.2.1 (single precision)
>>>> Reading file water.tpr, VERSION 3.2.1 (single precision)
>>>> How many groups do you want to calculate the RDF of?
>>>> 1
>>>> Select a reference group and 1 group
>>>> Group 0 ( OW) has 216 elements
>>>> There is one group in the index
>>>> There is one group in the index
>>>> trn version: GMX_trn_file (single precision)
>>>> Last frame 10 time 10.000
>>>> Segmentation fault
>>>> *******************************************************************************
>>>> 2. gromacs 3.1.5_pre2 on Mac OS X 10.3 (Panther). g_rdf gives the
>>>> following results (which appear to be the same as the above):
>>>> $ g_rdf -f water.trr -n oxygen.ndx -o rdf.xvg -s water.tpr
>>>> :-) G R O M A C S (-:
>>>> Go Rough, Oppose Many Angry Chinese Serial killers
>>>> :-) VERSION 3.1.5_pre2 (-:
>>>> (deleted some lines)
>>>> :-) g_rdf (-:
>>>> Option Filename Type Description
>>>> ------------------------------------------------------------
>>>> -f water.trr Input Generic trajectory: xtc trr trj
>>>> gro g96 pdb
>>>> -s water.tpr Input, Opt! Structure+mass(db): tpr tpb tpa
>>>> gro g96 pdb
>>>> -n oxygen.ndx Input, Opt! Index file
>>>> (deleted all remaining option lines)
>>>> Reading file water.tpr, VERSION 3.1.5_pre2 (single precision)
>>>> Reading file water.tpr, VERSION 3.1.5_pre2 (single precision)
>>>> How many groups do you want to calculate the RDF of?
>>>> 1
>>>> Select a reference group and 1 group
>>>> Group 0 ( OW) has 216 elements
>>>> There is one group in the index
>>>> There is one group in the index
>>>> trn version: GMX_trn_file
>>>> Last frame 10 time 10.000
>>>> Segmentation fault
>>>> *******************************************************************************
>>>> _______________________________________________
>>>> gmx-users mailing list gmx-users at gromacs.org
>>>> http://www.gromacs.org/mailman/listinfo/gmx-users
>>>> Please don't post (un)subscribe requests to the list. Use the www
>>>> interface or send it to gmx-users-request at gromacs.org.
>>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>>
>>>
>>>
>>> --David.
>>> ________________________________________________________________________
>>> David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
>>> Dept. of Cell and Molecular Biology, Uppsala University.
>>> Husargatan 3, Box 596, 75124 Uppsala, Sweden
>>> phone: 46 18 471 4205 fax: 46 18 511 755
>>> spoel at xray.bmc.uu.se spoel at gromacs.org http://folding.bmc.uu.se
>>> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>>> _______________________________________________
>>> gmx-users mailing list gmx-users at gromacs.org
>>> http://www.gromacs.org/mailman/listinfo/gmx-users
>>> Please don't post (un)subscribe requests to the list. Use the www
>>> interface or send it to gmx-users-request at gromacs.org.
>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>
>>
>>
>>
>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>> Dr. Kay-Eberhard Gottschalk
>> Applied Physics and Center for Nano Sciences
>> Ludwig-Maximilians University
>> Amalienstr. 54
>> 80799 Munich, Germany
>>
>> Phone: +49-89-2180 3436
>>
>>
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> gmx-users mailing list gmx-users at gromacs.org
>> http://www.gromacs.org/mailman/listinfo/gmx-users
>> Please don't post (un)subscribe requests to the list. Use the www
>> interface or send it to gmx-users-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
>
> _______________________________________________
> gmx-users mailing list gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
--
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596, 75124 Uppsala, Sweden
phone: 46 18 471 4205 fax: 46 18 511 755
spoel at xray.bmc.uu.se spoel at gromacs.org http://folding.bmc.uu.se
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
More information about the gromacs.org_gmx-users
mailing list