1#!/usr/local/bin/octave -qf
2#CATEGORY: Visualization and post-processing
3#SYNOPSIS: Create a PDB file from a JDFTx dry run
4
5function printUsage()
6	printf("\n\t Usage: dryRunToPDB <dryRunFile> [rep0=1] [rep1=1] [rep2=1]\n\n");
7	printf("\t Create <dryRunFile>.PDB and launch pymol to visualize the\n");
8	printf("\t system contained in the JDFTx dry-run output file <dryRunFile>.\n");
9	printf("\t Optionally specify repetitions along each lattice direction.\n");
10endfunction
11
12
13# Create PDB file from dry run output of JDFTx, with reps specifying repetitions along each of the axes
14function jdftxDryRunToPDB(outfile, reps)
15
16	#Get the lattice vectors:
17	command = sprintf("awk '$1==\"R\" && $2==\"=\" && Rdone!=1 {Rstart=NR; Rdone=1} Rstart && NR>Rstart && NR<=Rstart+3 { print $2,$3,$4}' %s", outfile);
18	fin = popen(command, "r");
19	R = reshape(fscanf(fin,"%f"), [3 3])';
20	pclose(fin);
21	Ra = R(:,1);
22	Rb = R(:,2);
23	Rc = R(:,3);
24
25	#Read the ion species names:
26	command = sprintf("awk '$1==\"ion\" {print $2}' %s", outfile);
27	fin = popen(command, "r");
28	nIons = 0;
29	while(!feof(fin))
30		nIons = nIons+1;
31		ionnames{nIons} = fgetl(fin);
32	endwhile
33	pclose(fin);
34
35	#Read the ion positions:
36	command = sprintf("awk '$1==\"ion\" {print $3, $4, $5}' %s", outfile);
37	fin = popen(command, "r");
38	ionpos = reshape(fscanf(fin,"%f"), 3, nIons);
39	pclose(fin);
40
41	#figure out coordinate system
42	command = sprintf("awk 'BEGIN {c=0} $1==\"coords-type\" && $2==\"cartesian\" {c++} END { print c}' %s", outfile);
43	fin = popen(command, "r");
44	cartesian = fscanf(fin,"%d");
45	pclose(fin);
46	if(cartesian>0)
47		ionpos = inv(R)*ionpos;
48	endif
49
50	#Output ions in PDB format, converted to cartesian coordinates, in Angstroms, with repetitions
51	fout = fopen(sprintf("%s.pdb",outfile), "w");
52	rnum=0;
53	for ia = 1:reps(1)
54	for ib = 1:reps(2)
55	for ic = 1:reps(3)
56	resName = sprintf("%1d%1d%1d",ia,ib,ic);
57	for i = 1:nIons
58		rnum = rnum + 1;
59		transformedpos = 0.5291772*((ia-1+ionpos(1,i))*Ra + (ib-1+ionpos(2,i))*Rb + (ic-1+ionpos(3,i))*Rc);
60		fprintf(fout, "HETATM%5d %4s %3s A%4d    %8.3f%8.3f%8.3f%6.2f%6.2f%10s%2s  \n", rnum, ionnames{i}, resName, rnum,
61				transformedpos(1), transformedpos(2), transformedpos(3), 0.0, 0.0, "", ionnames{i});
62	endfor
63	endfor
64	endfor
65	endfor
66	fclose(fout);
67
68	#For convenience, also invoke pymol to check the generated pdb:
69	fpml = fopen(sprintf("%s.pml",outfile), "w");
70	fprintf(fpml, "load %s.pdb\n", outfile);
71	fprintf(fpml, "show spheres, all\n");
72	fprintf(fpml, "set sphere_quality, 2\n");
73	fprintf(fpml, "set_view (\\\n");
74	fprintf(fpml, "   -0.877227128,   -0.234902740,    0.418676794,\\\n");
75	fprintf(fpml, "    0.474707872,   -0.554466426,    0.683535874,\\\n");
76	fprintf(fpml, "    0.071577206,    0.798368335,    0.597902119,\\\n");
77	fprintf(fpml, "    0.000024524,    0.000004083, -121.695564270,\\\n");
78	fprintf(fpml, "   12.041722298,   -0.376632214,    0.399373919,\\\n");
79	fprintf(fpml, "   91.902992249,  151.485549927,  -20.000000000 )\n");
80	fprintf(fpml, "center\n", outfile);
81	fprintf(fpml, "set ray_opaque_background, off\n", outfile);
82	fprintf(fpml, "ray 640,480\n", outfile);
83	fprintf(fpml, "png %s\n", outfile);
84	fclose(fpml);
85	system(sprintf("pymol %s.pml > /dev/null", outfile));
86endfunction
87
88arg_list = argv();
89if (nargin < 1 || strcmp(arg_list{1},"--help") || strcmp(arg_list{1},"-h"))
90	printUsage();
91	exit(1);
92endif
93
94dryRunFile = arg_list{1};
95reps = [1 1 1];
96if nargin >= 2; reps(1) = uint8(str2double(arg_list{2})); endif
97if nargin >= 3; reps(2) = uint8(str2double(arg_list{3})); endif
98if nargin >= 4; reps(3) = uint8(str2double(arg_list{4})); endif
99reps = max(reps, [1 1 1]);
100printf("Visualizing '%s' with %d x %d x %d repetitions.\n", dryRunFile, reps(1), reps(2), reps(3));
101jdftxDryRunToPDB(dryRunFile, reps);
102
103