!========================================================================= ! Program trg2atoms-03 : ! takes takes the atomic coordinates from a frame of a MD sequence ! and dumps all in an atom.inp ! it needs athe following files: ! - md.trg (ascii output trajectory from GULP ! - md.out output from GULP ! ! when lauched it requires 3 parameters: ! trg2atoms-03 a b c ! a = number of the frame to be processed ! b = edge to calculate 1=K, 2=L3 (no other are possible) ! c = Rmax to be used in the calculation ! ! can be compiled with g77 ! F. d'Acapito Ott 2011 ! !========================================================================= implicit real*8(a-h,o-z) character*15 flag character*38 cellpar character*100 BUFFER character achar*6, alphachar*12 character bchar*6, betachar*12 character cchar*6, gammachar*12 character latvec*41, flgnatoms*11 real mlatvec(3,3), milatvec(3,3), unity(3,3) character flgelem*38, dummy*80 character*2 elem(200), edge real x(200), y(200), z(200) real xo(200), yo(200), zo(200) real xf(200), yf(200), zf(200) ! input del numero del ciclo MD da trattare call GETARG(1,BUFFER) read(BUFFER,*) iframnum call GETARG(2,BUFFER) read(BUFFER,*) nedge call GETARG(3,BUFFER) read(BUFFER,*) rmax if (nedge.eq.1) edge='K ' if (nedge.eq.2) edge='L3' write(6,*) 'Program trg2atoms-03' write(6,*) iframnum open(3,file='md.trg') open(4,file='atoms.inp') open(5, file='md.out') open(7, file='atoms_org.inp') ! routine reading the cell parameters ! it uses the file md.out ! ! ! Getting the number of atoms in the cell do i=1,1000 read(5,'(a11)',END=10) flgnatoms if (flgnatoms.eq.' Formula =') then read(5,'(a45)') dummy read(5,'(a43, i3)') dummy, natoms goto 40 end if end do 40 write(6,*) natoms ! Reading the lattice vectors write(6,*) '# Cartesian lattice vectors Ang.' do i=1,1000 read(5,'(a41)',END=10) latvec if (latvec.eq.' Cartesian lattice vectors (Angstroms) :') then ! read(5,*) achar read (5,*) mlatvec(1,1),mlatvec(2,1),mlatvec(3,1) read (5,*) mlatvec(1,2),mlatvec(2,2),mlatvec(3,2) read (5,*) mlatvec(1,3),mlatvec(2,3),mlatvec(3,3) goto 30 end if end do 30 read(5,'(a38)',END=10) cellpar write(6,*) mlatvec(1,1), mlatvec(2,2),mlatvec(3,3) ! write the matrix as in the wikipedia formulas ! THE ELEments are indicated as aX where X runs from a to i! ! milatvec is the inverse matrix of mlatvec aa=mlatvec(1,1) ab=mlatvec(2,1) ac=mlatvec(3,1) ad=mlatvec(1,2) ae=mlatvec(2,2) af=mlatvec(3,2) ag=mlatvec(1,3) ah=mlatvec(2,3) ai=mlatvec(3,3) ! calculate the determinant deta = aa*ae*ai +ab*af*ag +ac*ad*ah -ac*ae*ag -af*ah*aa -ai*ab*ad ! calculate the inverse matrix milatvec(1,1)=(ae*ai -af*ah)/deta milatvec(1,2)=(af*ag -ad*ai)/deta milatvec(1,3)=(ad*ah -ae*ag)/deta milatvec(2,1)=(ac*ah -ab*ai)/deta milatvec(2,2)=(aa*ai -ac*ag)/deta milatvec(2,3)=(ab*ag -aa*ah)/deta milatvec(3,1)=(ab*af -ac*ae)/deta milatvec(3,2)=(ac*ad -aa*af)/deta milatvec(3,3)=(aa*ae -ab*ad)/deta write (6, *) 'Direct Matrix ' write(6,*) 'detA = ', deta write(6,'(f7.4, 2x,f7.4, 2x,f7.4)') mlatvec(1,1), mlatvec(2,1),mlatvec(3,1) write(6,'(f7.4, 2x,f7.4, 2x,f7.4)') mlatvec(1,2), mlatvec(2,2),mlatvec(3,2) write(6,'(f7.4, 2x,f7.4, 2x,f7.4)') mlatvec(1,3), mlatvec(2,3),mlatvec(3,3) write (6, *) 'Inverse Matrix ' write(6,'(f7.4, 2x,f7.4, 2x,f7.4)') milatvec(1,1), milatvec(2,1),milatvec(3,1) write(6,'(f7.4, 2x,f7.4, 2x,f7.4)') milatvec(1,2), milatvec(2,2),milatvec(3,2) write(6,'(f7.4, 2x,f7.4, 2x,f7.4)') milatvec(1,3), milatvec(2,3),milatvec(3,3) write(6, *) 'Check priduct matrix ' unity(1,1) = mlatvec(1,1)*milatvec(1,1) +mlatvec(2,1)*milatvec(1,2) +mlatvec(3,1)*milatvec(1,3) unity(2,1) = mlatvec(1,1)*milatvec(2,1) +mlatvec(2,1)*milatvec(2,2) +mlatvec(3,1)*milatvec(2,3) unity(3,1) = mlatvec(1,1)*milatvec(3,1) +mlatvec(2,1)*milatvec(3,2) +mlatvec(3,1)*milatvec(3,3) unity(1,2) = mlatvec(1,2)*milatvec(1,1) +mlatvec(2,2)*milatvec(1,2) +mlatvec(3,2)*milatvec(1,3) unity(2,2) = mlatvec(1,2)*milatvec(2,1) +mlatvec(2,2)*milatvec(2,2) +mlatvec(3,2)*milatvec(2,3) unity(3,2) = mlatvec(1,2)*milatvec(3,1) +mlatvec(2,2)*milatvec(3,2) +mlatvec(3,2)*milatvec(3,3) unity(1,3) = mlatvec(1,3)*milatvec(1,1) +mlatvec(2,3)*milatvec(1,2) +mlatvec(3,3)*milatvec(1,3) unity(2,3) = mlatvec(1,3)*milatvec(2,1) +mlatvec(2,3)*milatvec(2,2) +mlatvec(3,3)*milatvec(2,3) unity(3,3) = mlatvec(1,3)*milatvec(3,1) +mlatvec(2,3)*milatvec(3,2) +mlatvec(3,3)*milatvec(3,3) write(6, *) 'Product matrix' write(6,'(f6.4, 2x,f6.4, 2x,f6.4)') unity(1,1), unity(2,1),unity(3,1) write(6,'(f6.4, 2x,f6.4, 2x,f6.4)') unity(1,2), unity(2,2),unity(3,2) write(6,'(f6.4, 2x,f6.4, 2x,f6.4)') unity(1,3), unity(2,3),unity(3,3) ! REading the cell parameters write(6, *) 'Looking for the cell parameters' ! write(4,*) ' ' ! write(4, *) '# Cell parameters, a, b, c, alpha, beta, gamma' do i=1,1000 read(5,'(a38)',END=10) cellpar ! write(6,*) cellpar if (cellpar .eq. ' Cell parameters (Angstroms/Degrees):') then read(5, '(a6)') achar read (5, '(a6, F12.4, a12, F8.4)') achar, a, alphachar, alpha read (5, '(a6, F12.4, a12, F8.4)') bchar, b, betachar, beta read (5, '(a6, F12.4, a12, F8.4)') cchar, c, gammachar, gamma goto 20 end if end do 20 write(6, *) a, b, c, alpha, beta, gamma write (4,*) 'title = Frame num. ', iframnum write (4,'(a4, F7.4)') 'a = ', a write (4,'(a4, F7.4)') 'b = ', b write (4,'(a4, F7.4)') 'c = ', c write (7,*) 'title = Frame num. ', iframnum write (7,'(a4, F7.4)') 'a = ', a write (7,'(a4, F7.4)') 'b = ', b write (7,'(a4, F7.4)') 'c = ', c write (4,*) 'alpha = ', alpha write (4,*) 'beta = ', beta write (4,*) 'gamma = ', gamma write (4,*) 'Space = P1' write (7,*) 'alpha = ', alpha write (7,*) 'beta = ', beta write (7,*) 'gamma = ', gamma write (7,*) 'Space = P1' ! reading the chemical elements do i=1,1000 read(5,'(a38)',END=10) flgelem if (flgelem.eq.' Fractional coordinates of asymmetric') then do k=1,5 read(5, '(a45)') dummy end do do k=1,natoms read(5, '(a8, a2, a9, F8.6, a4, F8.6, a4, F8.6)') dummy, elem(k), dummy, & xo(k), dummy, yo(k), dummy, zo(k) ! write (6, *) dummy, elem(k) end do goto 50 end if end do 50 write (6, *) dummy write (4,*) 'core = ', elem(1) write(4,*) 'edge = ', edge write(4,*) 'rmax = ', rmax write (7,*) 'core = ', elem(1) write(7,*) 'edge = ', edge write(7,*) 'rmax = ', rmax ! routine reading the coordinates ! write (6, *) 'Cerco le coordinate...' ! write (4, *) a, b, c, alpha, beta, gamma write(4, *) ' ' write(4, *) 'atoms ' write(4, *) '! elem x y z tag & occ.' write(7, *) ' ' write(7, *) 'atoms ' write(7, *) '! elem x y z tag & occ.' read (3,*) version read (3,*) inat, ivar write(6,*) iframnum, inat k=0 do i=1,200000 read(3,'(a15)',END=10) flag ! write(6,*) flag if (flag.eq.'# Coordinates') then k=k+1 if (k.eq.iframnum) then do j=1,natoms read (3, *) x(j),y(j),z(j) ! calculate the fractional positions: dot procuct between coordinates and lattice vectors and then ! normalization xf(j) = ( x(j)*milatvec(1,1)+y(j)*milatvec(1,2)+z(j)*milatvec(1,3) ) if(xf(j).ge.1) xf(j)=xf(j)-1 if(xf(j).lt.0) xf(j)=xf(j)+1 yf(j) = ( x(j)*milatvec(2,1)+y(j)*milatvec(2,2)+z(j)*milatvec(2,3) ) if(yf(j).ge.1) yf(j)=yf(j)-1 if(yf(j).lt.0) yf(j)=yf(j)+1 zf(j) = ( x(j)*milatvec(3,1)+y(j)*milatvec(3,2)+z(j)*milatvec(3,3) ) if(zf(j).ge.1) zf(j)=zf(j)-1 if(zf(j).lt.0) zf(j)=zf(j)+1 write (4, '(a3, 3F12.8, a3, a2)') elem(j), xf(j),yf(j),zf(j), elem(j), ' 1' write (7, '(a3, 3F12.8, a3, a2)') elem(j), xo(j),yo(j),zo(j), elem(j), ' 1' ! write (4, *) elem(j), xf(j),yf(j),zf(j) end do goto 10 end if end if end do write(4,*) ' ' write(7,*) ' ' 10 close(3) close(4) close(5) close(7) end