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