1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20


21


22


23


24


25


26


27


28


29


30


31


32


33


34


35


36


37

#ifdef HAVE_CONFIG_H

38

#include <config.h>

39

#endif

40


41

#include <math.h>

42

#include <ctype.h>

43

#include "string2.h"

44

#include "sysstuff.h"

45

#include "typedefs.h"

46

#include "macros.h"

47

#include "vec.h"

48

#include "pbc.h"

49

#include "rmpbc.h"

50

#include "xvgr.h"

51

#include "copyrite.h"

52

#include "futil.h"

53

#include "statutil.h"

54

#include "tpxio.h"

55

#include "index.h"

56

#include "smalloc.h"

57

#include "fftgrid.h"

58

#include "calcgrid.h"

59

#include "nrnb.h"

60

#include "shift_util.h"

61

#include "pme.h"

62

#include "gstat.h"

63

#include "matio.h"

64


65

typedef struct

66

{

67

char *label;

68

int elem,mass;

69

real a[4], b[4], c;

70

} t_CM_table;

71


72


73


74


75


76


77


78

const t_CM_table CM_t[] =

79

{

80


81

{ "H", 1, 1, { 0.489918, 0.262003, 0.196767, 0.049879 },

82

{ 20.6593, 7.74039, 49.5519, 2.20159 },

83

0.001305 },

84

{ "C", 6, 12, { 2.26069, 1.56165, 1.05075, 0.839259 },

85

{ 22.6907, 0.656665, 9.75618, 55.5949 },

86

0.286977 },

87

{ "N", 7, 14, { 12.2126, 3.13220, 2.01250, 1.16630 },

88

{ 0.005700, 9.89330, 28.9975, 0.582600 },

89

11.529 },

90

{ "O", 8, 16, { 3.04850, 2.28680, 1.54630, 0.867000 },

91

{ 13.2771, 5.70110, 0.323900, 32.9089 },

92

0.250800 },

93

{ "Na", 11, 23, { 3.25650, 3.93620, 1.39980, 1.00320 },

94

{ 2.66710, 6.11530, 0.200100, 14.0390 },

95

0.404000 }

96

};

97


98

#define NCMT asize(CM_t)

99


100

typedef struct

101

{

102

int n_angles;

103

int n_groups;

104

double lambda;

105

double energy;

106

double momentum;

107

double ref_k;

108

double **F;

109

int nSteps;

110

int total_n_atoms;

111

} structure_factor;

112


113

typedef struct

114

{

115

rvec x;

116

int t;

117

} reduced_atom;

118


119

static void check_box_c(matrix box)

120

{

121

if (fabs(box[ZZ][XX]) > GMX_REAL_EPS*box[ZZ][ZZ] 

122

fabs(box[ZZ][YY]) > GMX_REAL_EPS*box[ZZ][ZZ])

123

gmx_fatal(FARGS,

124

"The last box vector is not parallel to the zaxis: %f %f %f",

125

box[ZZ][XX],box[ZZ][YY],box[ZZ][ZZ]);

126

}

127


128

static void do_rdf(char *fnNDX,char *fnTPS,char *fnTRX,

129

char *fnRDF,char *fnCNRDF, char *fnHQ,

130

bool bCM,bool bXY,bool bPBC,

131

real cutoff,real binwidth,real fade,int ng)

132

{

133

FILE *fp;

134

int status;

135

char outf1[STRLEN],outf2[STRLEN];

136

char title[STRLEN];

137

int g,natoms,i,j,k,nbin,j0,j1,n,nframes;

138

int **count;

139

char **grpname;

140

int *isize,isize_cm=0,nrdf=0,max_i;

141

atom_id **index,*index_cm=NULL;

142

#if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)

143

long long int *sum;

144

#else

145

double *sum;

146

#endif

147

real t,rmax2,cut2,r,r2,invbinw,normfac;

148

real segvol,spherevol,prev_spherevol,**rdf;

149

rvec *x,xcom,dx,*x_i1,xi;

150

real *inv_segvol,invvol,invvol_sum,rho;

151

bool *bExcl,bTop,bNonSelfExcl;

152

matrix box,box_pbc;

153

int **npairs,*nself;

154

atom_id ix,jx,***pairs;

155

t_topology top;

156

t_block *excl;

157

t_pbc pbc;

158


159

excl=NULL;

160


161

if (fnTPS) {

162

bTop=read_tps_conf(fnTPS,title,&top,&x,NULL,box,TRUE);

163

mk_single_top(&top);

164

if (bTop && !bCM)

165


166

excl=&(top.atoms.excl);

167

}

168

snew(grpname,ng+1);

169

snew(isize,ng+1);

170

snew(index,ng+1);

171

fprintf(stderr,"\nSelect a reference group and %d group%s\n",

172

ng,ng==1?"":"s");

173

if (fnTPS)

174

get_index(&top.atoms,fnNDX,ng+1,isize,index,grpname);

175

else

176

rd_index(fnNDX,ng+1,isize,index,grpname);

177


178

natoms=read_first_x(&status,fnTRX,&t,&x,box);

179

if ( !natoms )

180

gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");

181

if (fnTPS)

182


183

if ( natoms > top.atoms.nr )

184

gmx_fatal(FARGS,"Trajectory (%d atoms) does not match topology (%d atoms)",

185

natoms,top.atoms.nr);

186


187

for (i=0; i<=ng; i++)

188

for (j=0; j<isize[i]; j++)

189

if ( index[i][j] >= natoms )

190

gmx_fatal(FARGS,"Atom index (%d) in index group %s (%d atoms) larger "

191

"than number of atoms in trajectory (%d atoms)",

192

index[i][j],grpname[i],isize[i],natoms);

193


194

if (bCM) {

195


196

isize_cm=isize[0];

197

snew(index_cm,isize_cm);

198

for(i=0; i<isize[0]; i++)

199

index_cm[i]=index[0][i];

200

isize[0]=1;

201

index[0][0]=natoms;

202

srenew(index[0],isize[0]);

203


204

srenew(x,natoms+1);

205

}

206


207


208

copy_mat(box,box_pbc);

209

if (bXY) {

210

check_box_c(box);

211


212

box_pbc[ZZ][ZZ] = 2*max(box[XX][XX],box[YY][YY]);

213

}

214

rmax2 = 0.99*0.99*max_cutoff2(box_pbc);

215

nbin = (int)(sqrt(rmax2) / binwidth) + 1;

216

invbinw = 1.0 / binwidth;

217

cut2 = sqr(cutoff);

218


219

snew(count,ng);

220

snew(pairs,ng);

221

snew(npairs,ng);

222

snew(nself,ng);

223


224

snew(bExcl,natoms);

225

max_i = 0;

226

for(g=0; g<ng; g++) {

227

if (isize[g+1] > max_i)

228

max_i = isize[g+1];

229


230


231

snew(count[g],nbin+1);

232


233


234

snew(pairs[g],isize[0]);

235

snew(npairs[g],isize[0]);

236

for(i=0; i<isize[0]; i++) {

237

ix = index[0][i];

238

for(j=0; j < natoms; j++)

239

bExcl[j] = FALSE;

240


241

if (excl)

242

for( j = excl>index[ix]; j < excl>index[ix+1]; j++)

243

bExcl[excl>a[j]]=TRUE;

244

k = 0;

245

snew(pairs[g][i], isize[g+1]);

246

bNonSelfExcl = FALSE;

247

for(j=0; j<isize[g+1]; j++) {

248

jx = index[g+1][j];

249

if (!bExcl[jx])

250

pairs[g][i][k++]=jx;

251

else if (ix != jx)

252


253

bNonSelfExcl = TRUE;

254

if (ix == jx)

255

nself[g]++;

256

}

257

if (bNonSelfExcl) {

258

npairs[g][i]=k;

259

srenew(pairs[g][i],npairs[g][i]);

260

} else {

261


262

npairs[g][i]=1;

263

sfree(pairs[g][i]);

264

}

265

}

266

}

267

sfree(bExcl);

268


269

snew(x_i1,max_i);

270

nframes = 0;

271

invvol_sum = 0;

272

do {

273


274

copy_mat(box,box_pbc);

275

if (bPBC) {

276

rm_pbc(&top.idef,natoms,box,x,x);

277

if (bXY) {

278

check_box_c(box);

279

box_pbc[ZZ][ZZ] = 2*max(box[XX][XX],box[YY][YY]);

280

}

281

set_pbc(&pbc,box_pbc);

282


283

if (bXY)

284


285

box_pbc[ZZ][ZZ] = 1;

286

}

287

invvol = 1/det(box_pbc);

288

invvol_sum += invvol;

289


290

if (bCM) {

291


292

clear_rvec(xcom);

293

for(i=0; (i < isize_cm); i++) {

294

ix = index_cm[i];

295

rvec_inc(xcom,x[ix]);

296

}

297


298

for(j=0; (j<DIM); j++)

299

x[index[0][0]][j] = xcom[j] / isize_cm;

300

}

301


302

for(g=0; g<ng; g++) {

303


304

for(i=0; i<isize[g+1]; i++)

305

copy_rvec(x[index[g+1][i]],x_i1[i]);

306


307

for(i=0; i<isize[0]; i++) {

308

copy_rvec(x[index[0][i]],xi);

309

if (npairs[g][i] >= 0)

310


311

for(j=0; j<npairs[g][i]; j++) {

312

jx=pairs[g][i][j];

313

if (bPBC)

314

pbc_dx(&pbc,xi,x[jx],dx);

315

else

316

rvec_sub(xi,x[jx],dx);

317


318

if (bXY)

319

r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];

320

else

321

r2=iprod(dx,dx);

322

if (r2>cut2 && r2<=rmax2)

323

count[g][(int)(sqrt(r2)*invbinw)]++;

324

}

325

else

326


327

for(j=0; j<isize[g+1]; j++) {

328

if (bPBC)

329

pbc_dx(&pbc,xi,x_i1[j],dx);

330

else

331

rvec_sub(xi,x_i1[j],dx);

332

if (bXY)

333

r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];

334

else

335

r2=iprod(dx,dx);

336

if (r2>cut2 && r2<=rmax2)

337

count[g][(int)(sqrt(r2)*invbinw)]++;

338

}

339

}

340

}

341

nframes++;

342

} while (read_next_x(status,&t,natoms,x,box));

343

fprintf(stderr,"\n");

344


345

close_trj(status);

346


347

sfree(x);

348


349


350

invvol = invvol_sum/nframes;

351


352


353

snew(inv_segvol,nbin);

354

prev_spherevol=0;

355

for(i=0; (i<nbin); i++) {

356

r = (i+1)*binwidth;

357

if (bXY) {

358

spherevol=M_PI*r*r;

359

} else {

360

spherevol=(4.0/3.0)*M_PI*r*r*r;

361

}

362

segvol=spherevolprev_spherevol;

363

inv_segvol[i]=1.0/segvol;

364

prev_spherevol=spherevol;

365

}

366


367

snew(rdf,ng);

368

for(g=0; g<ng; g++) {

369


370

normfac = 1.0/(nframes*invvol*(isize[0]*isize[g+1]  nself[g]));

371


372


373

nrdf = max(nbin1,1+(2*fade/binwidth));

374

snew(rdf[g],nrdf);

375

for(i=0; (i<nbin1); i++) {

376

r = (i+0.5)*binwidth;

377

if ((fade > 0) && (r >= fade))

378

rdf[g][i] = 1+(count[g][i]*inv_segvol[i]*normfac1)*exp(16*sqr(r/fade1));

379

else

380

rdf[g][i] = count[g][i]*inv_segvol[i]*normfac;

381

}

382

for( ; (i<nrdf); i++)

383

rdf[g][i] = 1.0;

384

}

385


386

fp=xvgropen(fnRDF,"Radial Distribution","r","");

387

if (ng==1)

388

fprintf(fp,"@ subtitle \"%s%s\"\n",grpname[0],grpname[1]);

389

else {

390

fprintf(fp,"@ subtitle \"reference %s\"\n",grpname[0]);

391

xvgr_legend(fp,ng,grpname+1);

392

}

393

for(i=0; (i<nrdf); i++) {

394

fprintf(fp,"%10g", (i+0.5)*binwidth);

395

for(g=0; g<ng; g++)

396

fprintf(fp," %10g",rdf[g][i]);

397

fprintf(fp,"\n");

398

}

399

ffclose(fp);

400


401

do_view(fnRDF,NULL);

402


403


404

if (fnHQ) {

405

int nhq = 401;

406

real *hq,*integrand,Q;

407


408


409

rho = isize[1]*invvol;

410

snew(hq,nhq);

411

snew(integrand,nrdf);

412

for(i=0; (i<nhq); i++) {

413

Q = i*0.5;

414

integrand[0] = 0;

415

for(j=1; (j<nrdf); j++) {

416

r = (j+0.5)*binwidth;

417

integrand[j] = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);

418

integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]1.0);

419

}

420

hq[i] = print_and_integrate(debug,nrdf,binwidth,integrand,NULL,0);

421

}

422

fp=xvgropen(fnHQ,"h(Q)","Q(/nm)","h(Q)");

423

for(i=0; (i<nhq); i++)

424

fprintf(fp,"%10g %10g\n",i*0.5,hq[i]);

425

ffclose(fp);

426

do_view(fnHQ,NULL);

427

sfree(hq);

428

sfree(integrand);

429

}

430


431

if (fnCNRDF) {

432

normfac = 1.0/(isize[0]*nframes);

433

fp=xvgropen(fnCNRDF,"Cumulative Number RDF","r","number");

434

if (ng==1)

435

fprintf(fp,"@ subtitle \"%s%s\"\n",grpname[0],grpname[1]);

436

else {

437

fprintf(fp,"@ subtitle \"reference %s\"\n",grpname[0]);

438

xvgr_legend(fp,ng,grpname+1);

439

}

440

snew(sum,ng);

441

for(i=0; (i<nbin1); i++) {

442

fprintf(fp,"%10g",i*binwidth);

443

for(g=0; g<ng; g++) {

444

fprintf(fp," %10g",(real)((double)sum[g]*normfac));

445

sum[g] += count[g][i];

446

}

447

fprintf(fp,"\n");

448

}

449

ffclose(fp);

450

sfree(sum);

451


452

do_view(fnCNRDF,NULL);

453

}

454


455

for(g=0; g<ng; g++)

456

sfree(rdf[g]);

457

sfree(rdf);

458

}

459


460

typedef struct {

461

int ndata;

462

real kkk;

463

real intensity;

464

} t_xdata;

465


466

int comp_xdata(const void *a,const void *b)

467

{

468

t_xdata *xa,*xb;

469

real tmp;

470


471

xa = (t_xdata *)a;

472

xb = (t_xdata *)b;

473


474

if (xa>ndata == 0)

475

return 1;

476

else if (xb>ndata == 0)

477

return 1;

478

else {

479

tmp = xa>kkk  xb>kkk;

480

if (tmp < 0)

481

return 1;

482

else if (tmp > 0)

483

return 1;

484

else return 0;

485

}

486

}

487


488

static t_xdata *init_xdata(int nx,int ny)

489

{

490

int ix,iy,i,j,maxkx,maxky;

491

t_xdata *data;

492


493

maxkx = (nx+1)/2;

494

maxky = (ny+1)/2;

495

snew(data,nx*ny);

496

for(i=0; (i<nx); i++) {

497

for(j=0; (j<ny); j++) {

498

ix = abs((i < maxkx) ? i : (i  nx));

499

iy = abs((j < maxky) ? j : (j  ny));

500

data[ix*ny+iy].kkk = sqrt(ix*ix+iy*iy);

501

}

502

}

503

return data;

504

}

505


506

static void extract_sq(t_fftgrid *fftgrid,int nbin,real k_max,real lambda,

507

real count[],rvec box,int npixel,real *map[],

508

t_xdata data[])

509

{

510

int nx,ny,nz,nx2,ny2,nz2,la2,la12;

511

t_complex *ptr,*p0;

512

int i,j,k,maxkx,maxky,maxkz,n,ind,ix,iy;

513

real k1,kxy2,kz2,k2,z,kxy,kxy_max,cos_theta2,ttt,factor;

514

rvec lll,kk;

515


516


517


518


519

unpack_fftgrid(fftgrid,&nx,&ny,&nz,&nx2,&ny2,&nz2,

520

&la2,&la12,FALSE,(real **)&ptr);

521


522

maxkx = (nx+1)/2;

523

maxky = (ny+1)/2;

524

maxkz = nz/2+1;

525

factor = nbin/k_max;

526

for(i=0; (i<nx); i++) {

527

#define IDX(i,n) (i<=n/2) ? (i) : (in)

528

kk[XX] = IDX(i,nx);

529

for(j=0; (j<ny); j++) {

530

kk[YY] = IDX(j,ny);

531

ind = INDEX(i,j,0);

532

p0 = ptr + ind;

533

for(k=0; (k<maxkz); k++,p0++) {

534

if ((i==0) && (j==0) && (k==0))

535

continue;

536

kk[ZZ] = IDX(k,nz);

537

z = sqrt(sqr(p0>re)+sqr(p0>im));

538

kxy2 = sqr(kk[XX]) + sqr(kk[YY]);

539

k2 = kxy2+sqr(kk[ZZ]);

540

k1 = sqrt(k2);

541

ind = k1*factor;

542

if (ind < nbin) {

543


544


545


546


547


548


549


550


551


552


553

cos_theta2 = kxy2/k2;

554


555

ttt = z*0.5*(1+cos_theta2);

556

count[ind] += ttt;

557

ix = ((i < maxkx) ? i : (i  nx));

558

iy = ((j < maxky) ? j : (j  ny));

559

map[npixel/2+ix][npixel/2+iy] += ttt;

560

data[abs(ix)*ny+abs(iy)].ndata += 1;

561

data[abs(ix)*ny+abs(iy)].intensity += ttt;

562

}

563

else

564

fprintf(stderr,"k (%g) > k_max (%g)\n",k1,k_max);

565

}

566

}

567

}

568

}

569


570

typedef struct {

571

char *name;

572

int nelec;

573

} t_element;

574


575

static void do_sq(char *fnNDX,char *fnTPS,char *fnTRX,char *fnSQ,

576

char *fnXPM,real grid,real lambda,real distance,

577

int npixel,int nlevel)

578

{

579

FILE *fp;

580

t_element elem[] = { { "H", 1 }, { "C", 6 }, { "N", 7 }, { "O", 8 }, { "F", 9 }, { "S", 16 } };

581

#define NELEM asize(elem)

582

int status;

583

char title[STRLEN],*aname;

584

int natoms,i,j,k,nbin,j0,j1,n,nframes,pme_order;

585

real *count,**map;

586

char *grpname;

587

int isize;

588

atom_id *index;

589

real I0,C,t,k_max,factor,yfactor,segvol;

590

rvec *x,*xndx,box_size,kk,lll;

591

real fj0,*fj,max_spacing,r,lambda_1;

592

bool *bExcl;

593

matrix box;

594

int nx,ny,nz,nelectron;

595

atom_id ix,jx,**pairs;

596

splinevec *theta;

597

t_topology top;

598

t_fftgrid *fftgrid;

599

t_nrnb nrnb;

600

t_xdata *data;

601


602


603


604

fprintf(stderr,"\nSelect group for structure factor computation:\n");

605

get_index(&top.atoms,fnNDX,1,&isize,&index,&grpname);

606

natoms=read_first_x(&status,fnTRX,&t,&x,box);

607

if (isize < top.atoms.nr)

608

snew(xndx,isize);

609

else

610

xndx = x;

611

fprintf(stderr,"\n");

612


613

init_nrnb(&nrnb);

614


615

if ( !natoms )

616

gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");

617


618

if ( natoms > top.atoms.nr )

619

gmx_fatal(FARGS,"Trajectory (%d atoms) does not match topology (%d atoms)",

620

natoms,top.atoms.nr);

621


622


623

snew(fj,isize);

624

I0 = 0;

625

nelectron = 0;

626

for(i=0; (i<isize); i++) {

627

aname = *(top.atoms.atomname[index[i]]);

628

fj0 = 1;

629

if (top.atoms.atom[i].ptype == eptAtom) {

630

for(j=0; (j<NELEM); j++)

631

if (aname[0] == elem[j].name[0]) {

632

fj0 = elem[j].nelec;

633

break;

634

}

635

if (j == NELEM)

636

fprintf(stderr,"Warning: don't know number of electrons for atom %s\n",aname);

637

}

638


639

fj[i] = fj0  top.atoms.atom[index[i]].q;

640


641

nelectron += fj[i];

642


643

I0 += sqr(fj[i]);

644

}

645

if (debug) {

646


647

for(i=0; (i<isize); i++)

648

fprintf(debug,"Atom %3s%5d q = %10.4f f = %10.4f\n",

649

*(top.atoms.atomname[index[i]]),index[i],

650

top.atoms.atom[index[i]].q,fj[i]);

651

}

652


653


654

C = sqr(1.0/(ELECTRONMASS_keV*KILO*ELECTRONVOLT*1e7*distance));

655

fprintf(stderr,"C is %g\n",C);

656


657


658

nx = ny = nz = 0;

659

max_spacing = calc_grid(box,grid,&nx,&ny,&nz,1);

660

pme_order = max(4,1+(0.2/grid));

661

npixel = max(nx,ny);

662

data = init_xdata(nx,ny);

663


664

fprintf(stderr,"Largest grid spacing: %g nm, pme_order %d, %dx%d pixel on image\n",

665

max_spacing,pme_order,npixel,npixel);

666

fftgrid = init_pme(stdout,NULL,nx,ny,nz,pme_order,isize,FALSE,FALSE,eewg3D);

667


668


669

k_max = 1+sqrt(sqr(1+nx/2)+sqr(1+ny/2)+sqr(1+nz/2));

670


671


672

nbin = npixel;

673

snew(count,nbin+1);

674

snew(map,npixel);

675

for(i=0; (i<npixel); i++)

676

snew(map[i],npixel);

677


678

nframes = 0;

679

do {

680


681


682


683

for(j=0; (j<DIM); j++)

684

box_size[j] = box[j][j];

685


686


687

for(i=0; (i<isize); i++)

688

copy_rvec(x[index[i]],xndx[i]);

689


690


691

spread_on_grid(stdout,fftgrid,isize,pme_order,xndx,fj,box,FALSE,TRUE);

692


693


694

gmxfft3D(fftgrid,GMX_FFT_REAL_TO_COMPLEX,NULL);

695


696


697

extract_sq(fftgrid,nbin,k_max,lambda,count,box_size,npixel,map,data);

698


699

nframes++;

700

} while (read_next_x(status,&t,natoms,x,box));

701

fprintf(stderr,"\n");

702


703

close_trj(status);

704


705

sfree(x);

706


707


708

factor = k_max/(nbin);

709

yfactor = (1.0/nframes);

710

fp=xvgropen(fnSQ,"Structure Factor","q (1/nm)","S(q)");

711

fprintf(fp,"@ subtitle \"Lambda = %g nm. Grid spacing = %g nm\"\n",

712

lambda,grid);

713

factor *= lambda;

714

for(i=0; i<nbin; i++) {

715

r = (i+0.5)*factor*2*M_PI;

716

segvol = 4*M_PI*sqr(r)*factor;

717

fprintf(fp,"%10g %10g\n",r,count[i]*yfactor/segvol);

718

}

719

ffclose(fp);

720


721

do_view(fnSQ,NULL);

722


723

if (fnXPM) {

724

t_rgb rhi = { 0,0,0 }, rlo = { 1,1,1 };

725

real *tx,*ty,hi,inv_nframes;

726


727

hi = 0;

728

inv_nframes = 1.0/nframes;

729

snew(tx,npixel);

730

snew(ty,npixel);

731

for(i=0; (i<npixel); i++) {

732

tx[i] = inpixel/2;

733

ty[i] = inpixel/2;

734


735

for(j=0; (j<npixel); j++) {

736

map[i][j] *= inv_nframes;

737

hi = max(hi,map[i][j]);

738

}

739

}

740


741

fp = ffopen(fnXPM,"w");

742

write_xpm(fp,0,"Diffraction Image","Intensity","kx","ky",

743

nbin,nbin,tx,ty,map,0,hi,rlo,rhi,&nlevel);

744

fclose(fp);

745

sfree(tx);

746

sfree(ty);

747


748


749


750


751


752


753


754


755


756


757

}

758

}

759


760

t_complex *** rc_tensor_allocation(int x, int y, int z)

761

{

762

t_complex ***t;

763

int i,j;

764


765

snew(t,x);

766

t = (t_complex ***)calloc(x,sizeof(t_complex**));

767

if(!t) exit(fprintf(stderr,"\nallocation error"));

768

t[0] = (t_complex **)calloc(x*y,sizeof(t_complex*));

769

if(!t[0]) exit(fprintf(stderr,"\nallocation error"));

770

t[0][0] = (t_complex *)calloc(x*y*z,sizeof(t_complex));

771

if(!t[0][0]) exit(fprintf(stderr,"\nallocation error"));

772


773

for( j = 1 ; j < y ; j++)

774

t[0][j] = t[0][j1] + z;

775

for( i = 1 ; i < x ; i++) {

776

t[i] = t[i1] + y;

777

t[i][0] = t[i1][0] + y*z;

778

for( j = 1 ; j < y ; j++)

779

t[i][j] = t[i][j1] + z;

780

}

781

return t;

782

}

783


784

int return_atom_type (char *name)

785

{

786

typedef struct {

787

char *name;

788

int nh;

789

} t_united_h;

790

t_united_h uh[] = {

791

{ "CH1", 1 }, { "CH2", 2 }, { "CH3", 3 },

792

{ "CS1", 1 }, { "CS2", 2 }, { "CS3", 3 },

793

{ "CP1", 1 }, { "CP2", 2 }, { "CP3", 3 }

794

};

795

int i;

796


797

for(i0; (i<asize(uh)); i++)

798

if (strcmp(name,uh[i].name) == 0)

799

return NCMT1+uh[i].nh;

800


801

for(i=0; (i<NCMT); i++)

802

if (strncmp (name, CM_t[i].label,strlen(CM_t[i].label)) == 0)

803

return i;

804

gmx_fatal(FARGS,"\nError: atom (%s) not in list (%d types checked)!\n",

805

name,i);

806


807

return 0;

808

}

809


810

double CMSF (int type,int nh,double lambda, double sin_theta)

811


812


813


814


815


816

{

817

int i;

818

double tmp = 0.0, k2;

819


820


821


822


823


824

if (nh > 0) {

825

tmp = (CMSF (return_atom_type ("C"),0,lambda, sin_theta) +

826

nh*CMSF (return_atom_type ("H"),0,lambda, sin_theta));

827

}

828


829

else {

830

k2 = (sqr (sin_theta) / sqr (10.0 * lambda));

831

tmp = CM_t[type].c;

832

for (i = 0; (i < 4); i++)

833

tmp += CM_t[type].a[i] * exp (CM_t[type].b[i] * k2);

834

}

835

return tmp;

836

}

837


838

real **compute_scattering_factor_table (structure_factor * sf,int *nsftable)

839

{

840


841


842


843


844

int i, j;

845

double sin_theta,q,hc=1239.842;

846

real ** sf_table;

847


848


849

sf>momentum = ((double) (2. * 1000.0 * M_PI * sf>energy) / hc);

850

sf>lambda = hc / (1000.0 * sf>energy);

851

fprintf (stderr, "\nwavelenght = %f nm\n", sf>lambda);

852

*nsftable = NCMT+3;

853

snew (sf_table,*nsftable);

854

for (i = 0; (i < *nsftable); i++) {

855

snew (sf_table[i], sf>n_angles);

856

for (j = 0; j < sf>n_angles; j++) {

857

q = ((double) j * sf>ref_k);

858


859


860

sin_theta = q / (2.0 * sf>momentum);

861

if (i < NCMT)

862

sf_table[i][j] = CMSF (i,0,sf>lambda, sin_theta);

863

else

864

sf_table[i][j] = CMSF (i,iNCMT+1,sf>lambda, sin_theta);

865

}

866

}

867

return sf_table;

868

}

869


870

int * create_indexed_atom_type (reduced_atom * atm, int size)

871

{

872


873


874


875


876


877


878


879

int *index_atp, i, i_tmp, j;

880


881

snew (index_atp, 1);

882

i_tmp = 1;

883

index_atp[0] = atm[0].t;

884

for (i = 1; i < size; i++) {

885

for (j = 0; j < i_tmp; j++)

886

if (atm[i].t == index_atp[j])

887

break;

888

if (j == i_tmp) {

889

i_tmp++;

890

srenew (index_atp, i_tmp * sizeof (int));

891

index_atp[i_tmp  1] = atm[i].t;

892

}

893

}

894

i_tmp++;

895

srenew (index_atp, i_tmp * sizeof (int));

896

index_atp[i_tmp  1] = 0;

897

return index_atp;

898

}

899


900

void rearrange_atoms (reduced_atom * positions, t_trxframe *fr, atom_id * index,

901

int isize, t_topology * top, bool flag)

902


903

{

904

int i;

905


906

if (flag)

907

for (i = 0; i < isize; i++)

908

positions[i].t =

909

return_atom_type (*(top>atoms.atomname[index[i]]));

910

for (i = 0; i < isize; i++)

911

copy_rvec (fr>x[index[i]], positions[i].x);

912

}

913


914


915

int atp_size (int *index_atp)

916

{

917

int i = 0;

918


919

while (index_atp[i])

920

i++;

921

return i;

922

}

923


924

void compute_structure_factor (structure_factor * sf, matrix box,

925

reduced_atom * red, int isize, real start_q,

926

real end_q, int group,real **sf_table)

927

{

928

t_complex ***tmpSF;

929

rvec k_factor;

930

real kdotx, asf, kx, ky, kz, krr;

931

int kr, maxkx, maxky, maxkz, i, j, k, p, *counter;

932


933


934

k_factor[XX] = 2 * M_PI / box[XX][XX];

935

k_factor[YY] = 2 * M_PI / box[YY][YY];

936

k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ];

937


938

maxkx = (int) rint (end_q / k_factor[XX]);

939

maxky = (int) rint (end_q / k_factor[YY]);

940

maxkz = (int) rint (end_q / k_factor[ZZ]);

941


942

snew (counter, sf>n_angles);

943


944

tmpSF = rc_tensor_allocation(maxkx,maxky,maxkz);

945


946


947


948


949


950

fprintf(stderr,"\n");

951

for (i = 0; i < maxkx; i++) {

952

fprintf (stderr,"\rdone %3.1f%% ", (double)(100.0*(i+1))/maxkx);

953

kx = i * k_factor[XX];

954

for (j = 0; j < maxky; j++) {

955

ky = j * k_factor[YY];

956

for (k = 0; k < maxkz; k++)

957

if (i != 0  j != 0  k != 0) {

958

kz = k * k_factor[ZZ];

959

krr = sqrt (sqr (kx) + sqr (ky) + sqr (kz));

960

if (krr >= start_q && krr <= end_q) {

961

kr = (int) rint (krr/sf>ref_k);

962

if (kr < sf>n_angles) {

963

counter[kr]++;

964


965

for (p = 0; p < isize; p++) {

966


967

asf = sf_table[red[p].t][kr];

968


969

kdotx = kx * red[p].x[XX] +

970

ky * red[p].x[YY] + kz * red[p].x[ZZ];

971


972

tmpSF[i][j][k].re += cos (kdotx) * asf;

973

tmpSF[i][j][k].im += sin (kdotx) * asf;

974

}

975

}

976

}

977

}

978

}

979

}

980


981


982


983


984


985


986

for (i = 0; i < maxkx; i++) {

987

kx = i * k_factor[XX]; for (j = 0; j < maxky; j++) {

988

ky = j * k_factor[YY]; for (k = 0; k < maxkz; k++) {

989

kz = k * k_factor[ZZ]; krr = sqrt (sqr (kx) + sqr (ky)

990

+ sqr (kz)); if (krr >= start_q && krr <= end_q) {

991

kr = (int) rint (krr / sf>ref_k); if (kr <

992

sf>n_angles && counter[kr] != 0)

993

sf>F[group][kr] +=

994

(sqr (tmpSF[i][j][k].re) +

995

sqr (tmpSF[i][j][k].im))/ counter[kr];

996

}

997

}

998

}

999

} sfree (counter); free(tmpSF[0][0]); free(tmpSF[0]); free(tmpSF);

1000

}

1001


1002

void save_data (structure_factor * sf, char *file, int ngrps, real start_q,

1003

real end_q)

1004

{

1005


1006

FILE *fp;

1007

int i, g = 0;

1008

double *tmp, polarization_factor, A;

1009


1010

fp = xvgropen (file, "Scattering Intensity", "q (1/nm)",

1011

"Intensity (a.u.)");

1012


1013

snew (tmp, ngrps);

1014


1015

for (g = 0; g < ngrps; g++)

1016

for (i = 0; i < sf>n_angles; i++) {

1017


1018


1019


1020


1021


1022


1023


1024


1025

A = (double) (i * sf>ref_k) / (2.0 * sf>momentum);

1026

polarization_factor = 1  2.0 * sqr (A) * (1  sqr (A));

1027

sf>F[g][i] *= polarization_factor;

1028

}

1029

for (i = 0; i < sf>n_angles; i++) {

1030

if (i * sf>ref_k >= start_q && i * sf>ref_k <= end_q) {

1031

fprintf (fp, "%10.5f ", i * sf>ref_k);

1032

for (g = 0; g < ngrps; g++)

1033

fprintf (fp, " %10.5f ", (sf>F[g][i]) /( sf>total_n_atoms*

1034

sf>nSteps));

1035

fprintf (fp, "\n");

1036

}

1037

}

1038

ffclose (fp);

1039

}

1040


1041

int

1042

do_scattering_intensity (char* fnTPS, char* fnNDX, char* fnXVG, char *fnTRX,

1043

real start_q,real end_q, real energy,int ng)

1044

{

1045

int i,*isize,status,flags = TRX_READ_X,**index_atp;

1046

char **grpname,title[STRLEN];

1047

atom_id **index;

1048

t_topology top;

1049

t_trxframe fr;

1050

reduced_atom **red;

1051

structure_factor *sf;

1052

rvec *xtop;

1053

real **sf_table;

1054

int nsftable;

1055

matrix box;

1056

double r_tmp;

1057


1058

snew (sf, 1);

1059

sf>energy = energy;

1060


1061


1062

read_tps_conf (fnTPS, title, &top, &xtop, NULL, box, TRUE);

1063

sfree (xtop);

1064


1065


1066

snew (isize, ng);

1067

snew (index, ng);

1068

snew (grpname, ng);

1069


1070

fprintf (stderr, "\nSelect %d group%s\n", ng,

1071

ng == 1 ? "" : "s");

1072

if (fnTPS)

1073

get_index (&top.atoms, fnNDX, ng, isize, index, grpname);

1074

else

1075

rd_index (fnNDX, ng, isize, index, grpname);

1076


1077


1078

read_first_frame (&status, fnTRX, &fr, flags);

1079


1080

sf>total_n_atoms = fr.natoms;

1081


1082

snew (red, ng);

1083

snew (index_atp, ng);

1084


1085

r_tmp = max (box[XX][XX], box[YY][YY]);

1086

r_tmp = (double) max (box[ZZ][ZZ], r_tmp);

1087


1088

sf>ref_k = (2.0 * M_PI) / (r_tmp);

1089


1090

sf>n_angles = (int) rint (end_q / sf>ref_k);

1091


1092

snew (sf>F, ng);

1093

for (i = 0; i < ng; i++)

1094

snew (sf>F[i], sf>n_angles);

1095

for (i = 0; i < ng; i++) {

1096

snew (red[i], isize[i]);

1097

rearrange_atoms (red[i], &fr, index[i], isize[i], &top, TRUE);

1098

index_atp[i] = create_indexed_atom_type (red[i], isize[i]);

1099

}

1100

sf_table = compute_scattering_factor_table (sf,&nsftable);

1101


1102


1103

do {

1104

sf>nSteps++;

1105

for (i = 0; i < ng; i++) {

1106

rearrange_atoms (red[i], &fr, index[i], isize[i], &top,FALSE);

1107


1108

compute_structure_factor (sf, box, red[i], isize[i],

1109

start_q, end_q, i, sf_table);

1110

}

1111

}

1112

while (read_next_frame (status, &fr));

1113


1114

save_data (sf, fnXVG, ng, start_q, end_q);

1115


1116

return 0;

1117

}

1118


1119

int gmx_rdf(int argc,char *argv[])

1120

{

1121

static char *desc[] = {

1122

"The structure of liquids can be studied by either neutron or Xray",

1123

"scattering. The most common way to describe liquid structure is by a",

1124

"radial distribution function. However, this is not easy to obtain from",

1125

"a scattering experiment.[PAR]",

1126

"g_rdf calculates radial distribution functions in different ways.",

1127

"The normal method is around a (set of) particle(s), the other method",

1128

"is around the center of mass of a set of particles.",

1129

"With both methods rdf's can also be calculated around axes parallel",

1130

"to the zaxis with option [TT]xy[tt].[PAR]",

1131

"If a run input file is supplied ([TT]s[tt]), exclusions defined",

1132

"in that file are taken into account when calculating the rdf.",

1133

"The option [TT]cut[tt] is meant as an alternative way to avoid",

1134

"intramolecular peaks in the rdf plot.",

1135

"It is however better to supply a run input file with a higher number of",

1136

"exclusions. For eg. benzene a topology with nrexcl set to 5",

1137

"would eliminate all intramolecular contributions to the rdf.",

1138

"Note that all atoms in the selected groups are used, also the ones",

1139

"that don't have LennardJones interactions.[PAR]",

1140

"Option [TT]cn[tt] produces the cumulative number rdf.[PAR]"

1141

"To bridge the gap between theory and experiment structure factors can",

1142

"be computed (option [TT]sq[tt]). The algorithm uses FFT, the grid"

1143

"spacing of which is determined by option [TT]grid[tt]."

1144

};

1145

static bool bCM=FALSE,bXY=FALSE,bPBC=TRUE;

1146

static real cutoff=0,binwidth=0.002,grid=0.05,fade=0.0,lambda=0.1,distance=10;

1147

static int npixel=256,nlevel=20,ngroups=1;

1148

static real start_q=0.0, end_q=60.0, energy=12.0;

1149

t_pargs pa[] = {

1150

{ "bin", FALSE, etREAL, {&binwidth},

1151

"Binwidth (nm)" },

1152

{ "com", FALSE, etBOOL, {&bCM},

1153

"RDF with respect to the center of mass of first group" },

1154

{ "pbc", FALSE, etBOOL, {&bPBC},

1155

"Use periodic boundary conditions for computing distances" },

1156

{ "xy", FALSE, etBOOL, {&bXY},

1157

"Use only the x and y components of the distance" },

1158

{ "cut", FALSE, etREAL, {&cutoff},

1159

"Shortest distance (nm) to be considered"},

1160

{ "ng", FALSE, etINT, {&ngroups},

1161

"Number of secondary groups to compute RDFs around a central group" },

1162

{ "fade", FALSE, etREAL, {&fade},

1163

"From this distance onwards the RDF is tranformed by g'(r) = 1 + [g(r)1] exp((r/fade1)^2 to make it go to 1 smoothly. If fade is 0.0 nothing is done." },

1164

{ "grid", FALSE, etREAL, {&grid},

1165

"[HIDDEN]Grid spacing (in nm) for FFTs when computing structure factors" },

1166

{ "npixel", FALSE, etINT, {&npixel},

1167

"[HIDDEN]# pixels per edge of the square detector plate" },

1168

{ "nlevel", FALSE, etINT, {&nlevel},

1169

"Number of different colors in the diffraction image" },

1170

{ "distance", FALSE, etREAL, {&distance},

1171

"[HIDDEN]Distance (in cm) from the sample to the detector" },

1172

{ "wave", FALSE, etREAL, {&lambda},

1173

"[HIDDEN]Wavelength for Xrays/Neutrons for scattering. 0.1 nm corresponds to roughly 12 keV" },

1174


1175

{"startq", FALSE, etREAL, {&start_q},

1176

"Starting q (1/nm) "},

1177

{"endq", FALSE, etREAL, {&end_q},

1178

"Ending q (1/nm)"},

1179

{"energy", FALSE, etREAL, {&energy},

1180

"Energy of the incoming Xray (keV) "}

1181

};

1182

#define NPA asize(pa)

1183

char *fnTPS,*fnNDX;

1184

bool bSQ,bRDF;

1185


1186

t_filenm fnm[] = {

1187

{ efTRX, "f", NULL, ffREAD },

1188

{ efTPS, NULL, NULL, ffOPTRD },

1189

{ efNDX, NULL, NULL, ffOPTRD },

1190

{ efXVG, "o", "rdf", ffOPTWR },

1191

{ efXVG, "sq", "sq", ffOPTWR },

1192

{ efXVG, "cn", "rdf_cn", ffOPTWR },

1193

{ efXVG, "hq", "hq", ffOPTWR },

1194


1195

};

1196

#define NFILE asize(fnm)

1197


1198

CopyRight(stderr,argv[0]);

1199

parse_common_args(&argc,argv,PCA_CAN_VIEW  PCA_CAN_TIME  PCA_BE_NICE,

1200

NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);

1201


1202

fnTPS = ftp2fn_null(efTPS,NFILE,fnm);

1203

fnNDX = ftp2fn_null(efNDX,NFILE,fnm);

1204


1205

bSQ = opt2bSet("sq",NFILE,fnm);

1206

bRDF = opt2bSet("o",NFILE,fnm)  !bSQ;

1207


1208

if (bSQ) {

1209

if (!fnTPS)

1210

gmx_fatal(FARGS,"Need a tps file for calculating structure factors\n");

1211

}

1212

else {

1213

if (!fnTPS && !fnNDX)

1214

gmx_fatal(FARGS,"Neither index file nor topology file specified\n"

1215

" Nothing to do!");

1216

}

1217


1218

if (bSQ)

1219

do_scattering_intensity(fnTPS,fnNDX,opt2fn("sq",NFILE,fnm),ftp2fn(efTRX,NFILE,fnm),

1220

start_q, end_q, energy, ngroups );

1221


1222


1223


1224


1225

if (bRDF)

1226

do_rdf(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),

1227

opt2fn("o",NFILE,fnm),opt2fn_null("cn",NFILE,fnm),

1228

opt2fn_null("hq",NFILE,fnm),

1229

bCM,bXY,bPBC,cutoff,binwidth,fade,ngroups);

1230


1231

thanx(stderr);

1232


1233

return 0;

1234

}
