program geometry C *************************************************************** C Last modification: 10/10/2018 C Le os numeros atomicos na vertical C Last modification: 12/06/2018 (Paulo Pires) C Adicionada capacidade de leitura de arquivos *.dfr e ljname do DICE C Last modification: 18/03/2013 C Implementacao da leitura dos arquivos no formato do GROMACS C porem trajetoria com xyz C Last modification: 16/11/2011 C Implementacao da leitura dos arquivos no formato do GROMACS C Last modification: 08/02/2010 C O angulo de torcao foi modificado para ficar de 0 a 360. C Last modification: 14/12/2009 ERRADO o calculo da torcao C O angulo de torcao foi modificado para ficar de 0 a 360. C Last modification: 06/08/2008 C This program calculates distances, angles and torsions from C configurations of molecular systems. C Desenvolvido e implementado por Kaline Coutinho C *************************************************************** implicit none common / atomic / n, na, nan, m, natot, amass, aicon integer mnt, mnat, mnats, tabe parameter (mnt= 6, mnat= 7000, mnats= 60, tabe= 105) integer mnbond, mnangle, mndihed, mnctr parameter (mnbond= 100, mnangle= 300, mndihed= 300, mnctr= 6) real e, debye, rad2deg parameter (e= 18.2257, debye=0.2636, rad2deg=57.29577951) integer n(mnt), na(mnt), nan(mnat), m, natot real amass(tabe) character*2 aicon(tabe) integer iuin, iuxyz, iubnd, iuang, iudhd, iudip, iusol, iuall integer iuebnd, iueang, iuedhd, iueint, iueij, iulj integer iflag, mark1, mark2, mark3, mark4, mark5, mark6,ipt integer pot1, pot2, pot3, pot4, pot5, pot6 integer air(mnbond), ajr(mnbond), aic(mnctr) integer aia(mnangle),aja(mnangle),aka(mnangle) integer aid(mndihed),ajd(mndihed),akd(mndihed),ald(mndihed) integer nbnd, nang, ndhd, nimp, nchrg, nvdw, finput, ftinput integer ai, aj, ak, al, i, j, k, ni(3), nat, natomos, l1 integer ii, iunit, idec, icent, com, nconf, nctr, molec real rx(mnat), ry(mnat), rz(mnat) real qi(mnat), epsi(mnat), sigi(mnat), scl(mnats,mnats) real rij(mnbond), kr(mnbond), req(mnbond) real aijk(mnangle), kth(mnangle), theq(mnangle) real dijkl(mndihed), v(mndihed,3),f(mndihed,3) real eij(mnbond), eijk(mnangle), eijkl(mndihed) real dipx, dipy, dipz, dip, chrg, deltar, deltath, deltafi real distancia, angulo, torcao, teij, teijk, teijkl, te real lx, ly, lz, tchrg, lxm, lym, lzm, l5, l6 real fvdw14, fchg14, ftorsion, fvdw15, fchg15 real v12ij, v6ij, vcij, dij, qij, epsij, sigij real eb(mnats,mnats), eth(mnats,mnats), etr(mnats,mnats) real v12(mnats,mnats), v6(mnats,mnats), vc(mnats,mnats) real sr, sr2, sr6, sr12, sigma, epsilon, totij c real tmp, fi logical horario, einit, cargasdice, molecok character ty*2, rule*1 character cij(mnbond)*10, cijk(mnangle)*14, cijkl(mndihed)*18 character*1 cunit, cdec, ccent character*3 ci, cj, ck, cl character*35 namein, namexyz, nameout, ljname character*38 namesol, nameall character line*110, linei*82 iuin = 10 iuxyz = 11 iubnd = 12 iuang = 13 iudhd = 14 iudip = 15 iusol = 16 iuall = 17 iuebnd = 18 iueang = 19 iuedhd = 20 iueint = 21 iueij = 22 iulj = 23 fvdw14 = 1.00 fchg14 = 1.00 fvdw15 = 1.00 fchg15 = 1.00 ftorsion = 1.00 rule = '*' do i= 1, mnats do j= 1, mnats scl(i,j) = 1.00 enddo enddo write(6,'('' '')') write(6,*) " ****************************************************" write(6,*) " Program ANLGEOM " write(6,*) " Last modification: 12/06/2018 " write(6,*) " By Kaline Coutinho " write(6,*) "*****************************************************" write(6,'('' '')') write(6,'('' '')') write(6,'('' Keywords: '')') write(6,'(''format = 1=Gromacs 2=Tinker 3=Dice(default) '')') write(6,'(''nametop = Topology file of the solute '')') write(6,'('' Tinker=*.key Gromacs=*.top Dice=*.dfr'')') write(6,'(''nametraj = Trajectory file (*.xyz or *.pdb) '')') write(6,'(''ljname = Molecule informations (only for Dice)'')') write(6,'(''moltype = Number of types of molecules, m '')') write(6,'(''nmol = Quantity of each molecule type '')') write(6,'('' n(i),i=1,m '')') write(6,'(''natom = Number of atoms of each molecule type'')') write(6,'('' Atomic number of each molecule type '')') write(6,'('' in the same sequence that appears in '')') write(6,'('' trajectory file '')') write(6,'('' 8 1 1 (Example to water solvent) '')') c write(6,'(''internal = yes or no (To calculates E_internal) '')') write(6,'('' '')') write(6,'('' '')') c Le as keywords call input(namein, namexyz, finput, ftinput, einit, ljname) open(iuin,file=namein,status='old') write(6,*)'Openning: ', namein cargasdice = .false. natomos= 0 molec= 1 nctr = 0 nbnd = 0 nang = 0 ndhd = 0 nimp = 0 nconf= 0 nchrg= 0 nvdw = 0 iflag= 0 do while (iflag .eq. 0) read(iuin,'(a)',end=998) line if ((einit).and.(finput.eq.2)) then pot1 = index(line,'VDW-14-SCALE') pot2 = index(line,'CHG-14-SCALE') pot3 = index(line,'TORSIONUNIT') pot4 = index(line,'VDW-15-SCALE') pot5 = index(line,'CHG-15-SCALE') pot6 = index(line,'RADIUSRULE') if (pot1 .ne. 0) then pot1= pot1 + 12 read(line(pot1:),*) fvdw14 if (fvdw14 .gt. 1.00) fvdw14= 1.00/fvdw14 endif if (pot2 .ne. 0) then pot2= pot2 + 12 read(line(pot2:),*) fchg14 if (fchg14 .gt. 1.00) fchg14= 1.00/fchg14 endif if (pot3 .ne. 0) then pot3= pot3 + 12 read(line(pot3:),*) ftorsion if (ftorsion .gt. 1.00) ftorsion= 1.00/ftorsion endif if (pot4 .ne. 0) then pot4= pot4 + 12 read(line(pot4:),*) fvdw15 if (fvdw15 .gt. 1.00) fvdw15= 1.00/fvdw15 endif if (pot5 .ne. 0) then pot5= pot5 + 12 read(line(pot5:),*) fchg15 if (fchg15 .gt. 1.00) fchg15= 1.00/fchg15 endif if (pot6 .ne. 0) then if ((index(line,'GEOMETRIC')).eq. 0) rule='+' endif endif if (finput .eq. 2) then mark1= index(line,'BOND ') mark2= index(line,'ANGLE ') mark3= index(line,'TORSION ') mark4= index(line,'IMPTORS ') mark5= index(line,'CHARGE ') mark6= index(line,'VDW ') c le os atomos que estao ligados if (mark1 .ne. 0) then mark1= mark1 + 4 nbnd= nbnd + 1 if (nbnd .gt. mnbond) then write(6,*)'Error: Increase "mnbond" and recompiled' stop endif if (einit) then read(line(mark1:),*) ai, aj, kr(nbnd), req(nbnd) write(6,*) ' Read bond ', nbnd, ' between: ', ai, aj, : kr(nbnd),req(nbnd) else read(line(mark1:),*) ai, aj write(6,*) ' Read bond ', nbnd, ' between: ', ai, aj endif air(nbnd)= ai ajr(nbnd)= aj scl(ai,aj)= 0.00 scl(aj,ai)= 0.00 endif c le os atomos que formam angulos de ligacao if (mark2 .ne. 0) then mark2= mark2 + 5 nang= nang + 1 if (nang .gt. mnangle) then write(6,*)'Error: Increase "mnangle" and recompiled' stop endif if (einit) then read(line(mark2:),*) ai, aj, ak, kth(nang), theq(nang) write(6,*) ' Read angle ', nang,' between: ', ai, aj, ak, : kth(nang), theq(nang) else read(line(mark2:),*) ai, aj, ak write(6,*) ' Read angle ', nang,' between: ', ai, aj, ak endif aia(nang)= ai aja(nang)= aj aka(nang)= ak scl(ai,ak)= 0.00 scl(ak,ai)= 0.00 endif c le os atomos que formam angulos de torsao if (mark3 .ne. 0) then mark3= mark3 + 7 ndhd= ndhd + 1 if (ndhd .gt. mndihed) then write(6,*)'Error: Increase "mndihed" and recompiled' stop endif if (einit) then read(line(mark3:),*) ai, aj, ak, al, : (v(ndhd,i),f(ndhd,i),ni(i),i=1,3) write(6,*)' Read torsion ', ndhd,' between: ',ai,aj,ak, : al,(v(ndhd,i),f(ndhd,i),ni(i),i=1,3) do i= 1, 3 v(ndhd,i) = v(ndhd,i) * ftorsion enddo else read(line(mark3:),*) ai, aj, ak, al write(6,*) ' Read torsion ', ndhd, ' between: ',ai,aj,ak, : al endif aid(ndhd)= ai ajd(ndhd)= aj akd(ndhd)= ak ald(ndhd)= al scl(ai,al)= fvdw14 scl(al,ai)= fvdw14 c Busca o quinto vizinho do i= 1, nbnd if (air(i).eq.al) then if (scl(ai,ajr(i)).eq. 1.00) scl(ai,ajr(i))= fvdw15 elseif (ajr(i).eq.al) then if (scl(ai,air(i)).eq. 1.00) scl(ai,air(i))= fvdw15 elseif (air(i).eq.ai) then if (scl(al,ajr(i)).eq. 1.00) scl(al,ajr(i))= fvdw15 elseif (ajr(i).eq.ai) then if (scl(al,air(i)).eq. 1.00) scl(al,air(i))= fvdw15 endif enddo endif c le os atomos que formam angulos de torsao improprios if (mark4 .ne. 0) then mark4= mark4 + 7 nimp= nimp + 1 ndhd= ndhd + 1 if (ndhd .gt. mndihed) then write(6,*)'Error: Increase "mndihed" and recompiled' stop endif if (einit) then read(line(mark4:),*) ai, aj, ak, al, : (v(ndhd,i),f(ndhd,i),ni(i),i=1,3) write(6,*)' Read imptor ', nimp, ' between: ',ai,aj,ak, : al,(v(ndhd,i),f(ndhd,i),ni(i),i=1,3) else read(line(mark4:),*) ai, aj, ak, al write(6,*) ' Read imptor ', nimp, ' between: ',ai,aj,ak, : al endif aid(ndhd)= ai ajd(ndhd)= aj akd(ndhd)= ak ald(ndhd)= al endif c le cargas dos atomos if (mark5 .ne. 0) then mark5= mark5 + 6 nchrg= nchrg + 1 if (nchrg .gt. mnat) then write(6,*)'Error: Increase "mnat" and recompiled' stop endif read(line(mark5:),*) ai, chrg write(6,*) ' Read charge ', ai, ' : ', chrg qi(ai)= chrg * e endif c le parametros Lennard-Jones dos atomos if (mark6 .ne. 0) then mark6= mark6 + 4 nvdw= nvdw + 1 if (nvdw .gt. mnat) then write(6,*)'Error: Increase "mnat" and recompiled' stop endif read(line(mark6:),*) ai, sigma, epsilon write(6,*) ' Read Lennard-Jones ', ai,' : ',sigma,epsilon epsi(ai)= 2.0 * sqrt(epsilon) if (rule .eq. '*') then sigi(ai)= sqrt(sigma) else sigi(ai)= sigma/2.0 endif endif c leitura de arquivo *.top do gromacs linhas 323 a 499 elseif (finput .eq. 1) then mark1= index(line,'[ atoms ]') if (mark1 .ne. 0) then tchrg= 0.00 read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,'(i6,39x,f11.3)',err=100,end=100) ai, chrg if (ai .gt. 0) then tchrg = tchrg + chrg nchrg= nchrg + 1 if (nchrg .gt. mnat) then write(6,*)'Error: Increase "mnat" and recompiled' stop endif qi(ai)= chrg * e i = 1 write(6,'(a,i6,3a,f11.3)') ' Read charge ',ai,' ', : aicon(nan(ai)), ' : ', chrg read(iuin,'(a)',end=998) line else write(6,*) 'Total charge: ', tchrg goto 100 endif enddo 100 continue endif mark1= index(line,'[ center ]') if (mark1 .ne. 0) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=110,end=110) nctr if (nctr .gt. mnctr) then write(6,*)'Error: Increase "mnctr" and recompiled' stop endif read(iuin,*) (aic(j),j=1,nctr) write(6,*) ' Read center ', nctr, ' between: ', : (aic(j),j=1,nctr) read(iuin,'(a)',end=998) line enddo 110 continue endif mark2= index(line,'[ bonds ]') if (mark2 .ne. 0) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=200,end=200) ai, aj nbnd= nbnd + 1 write(6,*) ' Read bond ', nbnd, ' between: ', ai, aj if (nbnd .gt. mnbond) then write(6,*)'Error: Increase "mnbond" and recompiled' stop endif air(nbnd)= ai ajr(nbnd)= aj read(iuin,'(a)',end=998) line enddo 200 continue endif mark3= index(line,'[ angles ]') if (mark3 .ne. 0) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=300,end=300) ai, aj, ak nang= nang + 1 write(6,*) ' Read angle ', nang,' between: ', ai, aj, ak if (nang .gt. mnangle) then write(6,*)'Error: Increase "mnangle" and recompiled' stop endif aia(nang)= ai aja(nang)= aj aka(nang)= ak read(iuin,'(a)',end=998) line enddo 300 continue endif mark4= index(line,'[ dihedrals ]') if ((mark4 .ne. 0).and.(ndhd .eq. 0)) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=400,end=400) ai, aj, ak, al ndhd= ndhd + 1 write(6,*) ' Read torsion ',ndhd,' between: ',ai,aj,ak,al if (ndhd .gt. mndihed) then write(6,*)'Error: Increase "mndihed" and recompiled' stop endif aid(ndhd)= ai ajd(ndhd)= aj akd(ndhd)= ak ald(ndhd)= al read(iuin,'(a)',end=998) line enddo 400 continue endif mark5= index(line,'[ dihedrals ]') if ((mark5 .ne. 0).and.(ndhd .gt. 0)) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=500,end=500) ai, aj, ak, al nimp= nimp + 1 ndhd= ndhd + 1 write(6,*) ' Read imptor ', nimp,' between: ',ai,aj,ak,al if (ndhd .gt. mndihed) then write(6,*)'Error: Increase "mndihed" and recompiled' stop endif aid(ndhd)= ai ajd(ndhd)= aj akd(ndhd)= ak ald(ndhd)= al read(iuin,'(a)',end=998) line enddo 500 continue endif elseif (finput .eq. 3) then c leitura de arquivo ljname do Dice if (.not. cargasdice) then open(iulj,file=ljname,status='old') write(6,*)'Openning: ', ljname read(iulj,*) read(iulj,'(I2)') l1 write(6,*)'File ',ljname,' has ',l1,' molecules' molecok = .false. do while (.not. molecok) molecok = .true. read(iulj,'(I2)',err=504,end=504) natomos if (natomos .ne. na(molec)) then do 503 i = 1, natomos read(iulj,*) 503 continue molec = molec + 1 molecok = .false. endif enddo 504 if (.not. molecok) then write(6,*)'Error: molecules in ', ljname,' have different : number of atoms than those defined in the input file' stop endif write(6,*)'Using ONLY data of molecule ',molec,'.' tchrg= 0.00 do 505 i = 1, natomos read(iulj,*,err=505,end=505) l1, ai, lx, ly, lz, chrg, l5, : l6 if (ai .gt. 0) then tchrg = tchrg + chrg nchrg= nchrg + 1 if (nchrg .gt. mnat) then write(6,*)'Error: Increase "mnat" and recompiled' stop endif qi(i)= chrg * e write(6,'(a,i6,3a,f11.3)') ' Read charge ',i,' ', : aicon(nan(i)), ' : ', chrg endif 505 continue close(iulj) write(6,*) 'File ',ljname,' closed.' write(6,*) 'Total charge: ', tchrg cargasdice= .true. endif c leitura de arquivo *.fdr do dice: linhas 501 a 615 mark2= index(line,'$bond') if (mark2 .ne. 0) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=600,end=600) ai, aj nbnd= nbnd + 1 write(6,*) ' Read bond ', nbnd, ' between: ', ai, aj if (nbnd .gt. mnbond) then write(6,*)'Error: Increase "mnbond" and recompiled' stop endif air(nbnd)= ai ajr(nbnd)= aj read(iuin,'(a)',end=998) line enddo 600 continue endif mark3= index(line,'$angle') if (mark3 .ne. 0) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=610,end=610) ai, aj, ak nang= nang + 1 write(6,*) ' Read angle ', nang,' between: ', ai, aj, ak if (nang .gt. mnangle) then write(6,*)'Error: Increase "mnangle" and recompiled' stop endif aia(nang)= ai aja(nang)= aj aka(nang)= ak read(iuin,'(a)',end=998) line enddo 610 continue endif mark4= index(line,'$dihedral') if ((mark4 .ne. 0).and.(ndhd .eq. 0)) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=620,end=620) ai, aj, ak, al ndhd= ndhd + 1 write(6,*) ' Read torsion ',ndhd,' between: ',ai,aj,ak,al if (ndhd .gt. mndihed) then write(6,*)'Error: Increase "mndihed" and recompiled' stop endif aid(ndhd)= ai ajd(ndhd)= aj akd(ndhd)= ak ald(ndhd)= al read(iuin,'(a)',end=998) line enddo 620 continue endif mark5= index(line,'$improper dihedral') if ((mark5 .ne. 0).and.(ndhd .gt. 0)) then read(iuin,'(a)',end=998) line com = 0 if (line(1:1).eq.';') com = 1 do while (com .eq. 1) read(iuin,'(a)',end=998) line if (line(1:1).eq.';') then com = 1 else com = 0 endif enddo do while (com .eq. 0) read(line,*,err=630,end=630) ai, aj, ak, al nimp= nimp + 1 ndhd= ndhd + 1 write(6,*) ' Read imptor ', nimp,' between: ',ai,aj,ak,al if (ndhd .gt. mndihed) then write(6,*)'Error: Increase "mndihed" and recompiled' stop endif aid(ndhd)= ai ajd(ndhd)= aj akd(ndhd)= ak ald(ndhd)= al read(iuin,'(a)',end=998) line enddo 630 continue endif endif enddo 998 continue close(iuin) write(6,*)'Closing: ', namein write(6,*) if ((einit).and.(finput.eq.2)) then write(6,*) 'Intramolecular scale factors: ' write(6,*) ' EPSILONRULE = GEOMETRIC' if (rule .eq. '*') then write(6,*) ' RADIUSRULE = GEOMETRIC' else write(6,*) ' RADIUSRULE = ARITMETRIC' endif write(6,*) ' VDW-14-SCALE = ', fvdw14 write(6,*) ' CHG-14-SCALE = ', fvdw14 write(6,*) ' TORSIONUNIT = ', ftorsion write(6,*) ' VDW-15-SCALE = ', fvdw15 write(6,*) ' CHG-15-SCALE = ', fvdw15 write(6,*) endif c write(6,*) 'Connectivity between atoms: ' c do i= 1, na(1)-1 c write(6,*) i,': ', (j,scl(i,j),j=i+1,na(1)) c enddo write(6,*) write(6,*) 'Summary of topology: ' write(6,*) nbnd,' bond distances will be analyzed.' write(6,*) nang,' bond angles will be analyzed.' write(6,*) nimp,' improper torsion angles will be analyzed.' write(6,*) ndhd-nimp,' torsion angles will be analyzed.' write(6,*) nchrg,' read atom charges.' write(6,*) c gera cabecalho do arquivo *.bnd do i=1,nbnd do j= 1, 2 if (j.eq.1) ii= air(i) if (j.eq.2) ii= ajr(i) if (ii.le. 9) then icent= 0 idec = 0 iunit= ii elseif (ii .le. 99) then icent= 0 idec = int(ii/10) iunit= ii - idec*10 elseif (ii .le. 999) then icent = int(ii/100) idec = int((ii-icent*100)/10) iunit= ii - (icent*100 + idec*10) endif ccent = char(icent+48) cdec = char(idec +48) cunit = char(iunit+48) if (j.eq.1) ci= ccent // cdec // cunit if (j.eq.2) cj= ccent // cdec // cunit enddo cij(i)='(' // ci // '-' // cj // ') ' enddo c gera cabecalho do arquivo *.ang do i=1,nang do j= 1, 3 if (j.eq.1) ii= aia(i) if (j.eq.2) ii= aja(i) if (j.eq.3) ii= aka(i) if (ii.le. 9) then icent= 0 idec = 0 iunit= ii elseif (ii .le. 99) then icent= 0 idec = int(ii/10) iunit= ii - idec*10 elseif (ii .le. 999) then icent = int(ii/100) idec = int((ii-icent*100)/10) iunit= ii - (icent*100 + idec*10) endif ccent = char(icent+48) cdec = char(idec +48) cunit = char(iunit+48) if (j.eq.1) ci= ccent // cdec // cunit if (j.eq.2) cj= ccent // cdec // cunit if (j.eq.3) ck= ccent // cdec // cunit enddo cijk(i)='(' // ci // '-' // cj // '-' // ck //') ' enddo c gera cabecalho do arquivo *.dhd do i=1,ndhd do j= 1, 4 if (j.eq.1) ii= aid(i) if (j.eq.2) ii= ajd(i) if (j.eq.3) ii= akd(i) if (j.eq.4) ii= ald(i) if (ii.le. 9) then icent= 0 idec = 0 iunit= ii elseif (ii .le. 99) then icent= 0 idec = int(ii/10) iunit= ii - idec*10 elseif (ii .le. 999) then icent = int(ii/100) idec = int((ii-icent*100)/10) iunit= ii - (icent*100 + idec*10) endif ccent = char(icent+48) cdec = char(idec +48) cunit = char(iunit+48) if (j.eq.1) ci= ccent // cdec // cunit if (j.eq.2) cj= ccent // cdec // cunit if (j.eq.3) ck= ccent // cdec // cunit if (j.eq.4) cl= ccent // cdec // cunit enddo cijkl(i)='('// ci // '-' // cj // '-' // ck // '-' // cl //') ' enddo c gerar o nome do arquivo de saida ipt= index(namexyz,'.') ipt= ipt - 1 c escreve nos arquivo *.bnd nameout= namexyz(:ipt) // '-bnd.avr' open(iubnd,file=nameout, status='unknown') write(iubnd,*)' NMOVE ',(cij(i),i=1,nbnd) c escreve nos arquivo *.ang nameout= namexyz(:ipt) // '-ang.avr' open(iuang,file=nameout, status='unknown') write(iuang,*)' NMOVE ',(cijk(i),i=1,nang) c escreve nos arquivo *.dhd nameout= namexyz(:ipt) // '-dhd.avr' open(iudhd,file=nameout, status='unknown') write(iudhd,*)' NMOVE ', (cijkl(i),i=1,ndhd) c escreve nos arquivo *.dip if format .ne. dice (Paulo Pires) nameout= namexyz(:ipt) // '-dip.avr' open(iudip,file=nameout, status='unknown') write(iudip,*)' NMOVE Dip_x Dip_y Dip_z Dipole_tot' c escreve nos arquivos de energia *.ebnd, *.eang, *.edhb, *.eint, c *.eijall if (einit) then nameout= namexyz(:ipt) // 'ebnd' open(iuebnd,file=nameout, status='unknown') write(iuebnd,*)' NMOVE ',(cij(i),i=1,nbnd) nameout= namexyz(:ipt) // 'eang' open(iueang,file=nameout, status='unknown') write(iueang,*)' NMOVE ', (cijk(i),i=1,nang) nameout= namexyz(:ipt) // 'edhd' open(iuedhd,file=nameout, status='unknown') write(iuedhd,*)' NMOVE ', (cijkl(i),i=1,ndhd) nameout= namexyz(:ipt) // 'eint' open(iueint,file=nameout, status='unknown') write(iueint,*)' NMOVE Eij Eijk Eijkl E12 -E6', : ' Eq Etotal' nameout= namexyz(:ipt) // 'eijall' open(iueij,file=nameout, status='unknown') write(iueij,*) ' NMOVE i j Ebnd Eth Etor ', : 'E12 -E6 Eq Etotal' endif c escreve nos arquivo _sol.xyz ipt=ipt-1 namesol= namexyz(:ipt) // '_sol.xyz' open(iusol,file=namesol, status= 'unknown') if (finput .ne. 3) then nameall= namexyz(:ipt) // '_all.xyz' open(iuall,file=nameall, status= 'unknown') endif open(iuxyz,file=namexyz,status='old') write(6,*)'Openning: ', namexyz nconf= 0 do while (iflag .eq. 0) c loop sobre todas as configuracoes dipx = 0.0 dipy = 0.0 dipz = 0.0 nconf= nconf + 1 if (ftinput .eq. 2) then c le o format xyx geral read(iuxyz,*,end=9999) nat if (nat .ne. natot) then write(6,*) 'Error: number of atoms in xyz file is ', : 'different of "natot" keyword' stop endif write(iusol,'(i8)') na(1) if (finput .ne. 3) write(iuall,'(i8)') natot read(iuxyz,'(a)') linei mark1 = index(linei,'L =') if (mark1 .ne.0) then mark1= mark1 + 4 read(linei(mark1:),*) lx, ly, lz write(iusol,'(a,i7,a,3f9.4)') : ' Configuration number : ', nconf, 'L = ', lx, ly, lz if (finput .ne. 3) write(iuall,'(a,i7,a,3f9.4)') : ' Configuration number : ', nconf, 'L = ', lx, ly, lz else if (finput .ne. 3) write(iuall,'('' Configuration number : : '',i7)') nconf endif do i= 1, nat read(iuxyz,'(a)')linei j=1 k= j+ 1 do while(linei(j:k).eq.' ') j=j+1 k=k+1 enddo j=j+1 k=k+1 read(linei(j:k),*) ty do while(linei(j:k).eq.' ') j=j+1 k=k+1 enddo j=j+4 read(linei(j:),*) rx(i), ry(i), rz(i) if (i .le. na(1)) : write(iusol,'(2x,a,3(x,f11.5))') ty, rx(i), ry(i), rz(i) if (finput .ne. 3) write(iuall,'(2x,a,3(x,f11.5))') ty, : rx(i), ry(i), rz(i) enddo c calcula o centro dos atomos do center if (nctr .gt. 0) then k= na(1)+1 rx(k)= 0.00 ry(k)= 0.00 rz(k)= 0.00 do j= 1, nctr write(6,*) j, rx(aic(j)), ry(aic(j)), rz(aic(j)) rx(k) = rx(k) + rx(aic(j)) ry(k) = ry(k) + ry(aic(j)) rz(k) = rz(k) + rz(aic(j)) enddo rx(k) = rx(k)/real(nctr) ry(k) = ry(k)/real(nctr) rz(k) = rz(k)/real(nctr) write(6,*) k, rx(k), ry(k), rz(k) k = 28 write(6,*) k, rx(k), ry(k), rz(k) endif else c le o format pdb mark1 = 0 read(iuxyz,'(a)',end=9999)linei write(iusol,'(i8)') na(1) if (finput .ne. 3) write(iuall,'(i8)') natot mark1 = index(linei,'CRYST1') do while (mark1 .eq. 0) read(iuxyz,'(a)')linei mark1 = index(linei,'CRYST1') enddo read(linei(7:),*) lx, ly, lz lxm = lx / 2.0 lym = ly / 2.0 lzm = lz / 2.0 write(iusol,'(a,i7,a,3f9.4)') : ' Configuration number : ', nconf, 'L = ', lx, ly, lz if (finput .ne. 3) write(iuall,'(a,i7,a,3f9.4)') : ' Configuration number : ', nconf, 'L = ', lx, ly, lz read(iuxyz,'(a)')linei do i= 1, natot read(iuxyz,'(30x,3f8.3)') rx(i), ry(i), rz(i) rx(i) = rx(i) - lxm ry(i) = ry(i) - lym rz(i) = rz(i) - lzm if (i .le. na(1)) : write(iusol,'(2x,a,3(x,f11.5))') aicon(nan(i)), rx(i), : ry(i), rz(i) if (finput .ne. 3) write(iuall,'(2x,a,3(x,f11.5))') : aicon(nan(i)), rx(i), ry(i), rz(i) enddo read(iuxyz,'(a)')linei read(iuxyz,'(a)')linei endif write(6,*)' Read configuration: ', nconf c calculando as distancias e energias de ligacao teij = 0.00 deltar= 0.00 do i= 1, nbnd rij(i)= distancia(rx(air(i)),ry(air(i)),rz(air(i)), : rx(ajr(i)),ry(ajr(i)),rz(ajr(i))) if ((einit).and.(finput.eq.2)) then deltar = req(i)-rij(i) eij(i) = kr(i)*deltar*deltar teij = teij + eij(i) eb(air(i),ajr(i))= eij(i) eb(ajr(i),air(i))= eij(i) endif enddo write(iubnd,*) nconf,(rij(i),i=1,nbnd) if (einit) write(iuebnd,*) nconf,(eij(i),i=1,nbnd) c calculando os angulos e energias de ligacao teijk = 0.00 deltath= 0.00 do i= 1, nang aijk(i)= angulo(rx(aia(i)),ry(aia(i)),rz(aia(i)), : rx(aja(i)),ry(aja(i)),rz(aja(i)), : rx(aka(i)),ry(aka(i)),rz(aka(i))) if ((einit).and.(finput.eq.2)) then deltath= (theq(i)-aijk(i))/rad2deg eijk(i)= kth(i)*deltath*deltath teijk= teijk + eijk(i) eth(aia(i),aka(i))= eijk(i) eth(aka(i),aia(i))= eijk(i) endif enddo write(iuang,*) nconf,(aijk(i),i=1,nang) if (einit) write(iueang,*) nconf,(eijk(i),i=1,nang) c calculando os angulos e energias de torsao teijkl = 0.00 do i= 1, ndhd dijkl(i)= torcao(rx(aid(i)),ry(aid(i)),rz(aid(i)), : rx(ajd(i)),ry(ajd(i)),rz(ajd(i)), : rx(akd(i)),ry(akd(i)),rz(akd(i)), : rx(ald(i)),ry(ald(i)),rz(ald(i))) c determina o sentido da rotacao torcional if (horario(rx(aid(i)),ry(aid(i)),rz(aid(i)), : rx(ajd(i)),ry(ajd(i)),rz(ajd(i)), : rx(akd(i)),ry(akd(i)),rz(akd(i)), : rx(ald(i)),ry(ald(i)),rz(ald(i)))) : dijkl(i)= -dijkl(i) if ((einit).and.(finput.eq.2)) then eijkl(i)= 0.00 c tmp= 0.00 do j= 1, 3 deltafi= (j*dijkl(i)+f(i,j))/rad2deg eijkl(i)= eijkl(i)+ v(i,j)*(1.0+cos(deltafi)) teijkl= teijkl + eijkl(i) enddo c fi= dijkl(i)/rad2deg c tmp= v(i,1)*(1.0+cos(fi)) + c : v(i,2)*(1.0-cos(2*fi)) + c : v(i,3)*(1.0+cos(3*fi)) c if ((eijkl(i)-tmp).gt.1E-04) then c write(6,*) 'Problem in the torsion energy:',eijkl(i), tmp c stop c endif etr(aid(i),ald(i))= eijkl(i) etr(ald(i),aid(i))= eijkl(i) endif enddo write(iudhd,*) nconf,(dijkl(i),i=1,ndhd) if (einit) write(iuedhd,*) nconf,(eijkl(i),i=1,ndhd) c calculo da energia Lennard-Jones e Coulomb if ((einit).and.(finput.eq.2)) then v12ij = 0.00 v6ij = 0.00 vcij = 0.00 do i= 1, na(1)-1 do j= i+1, na(1) if (scl(i,j) .gt. 0.00) then dij = distancia(rx(i),ry(i),rz(i),rx(j),ry(j),rz(j)) qij = qi(i) * qi(j) epsij = epsi(i) * epsi(j) if (rule .eq. '*') then sigij = sigi(i) * sigi(j) else sigij = sigi(i) + sigi(j) endif sr = sigij / dij sr2 = sr * sr sr6 = sr2 * sr2 * sr2 sr12 = sr6 * sr6 v12(i,j)= scl(i,j)*epsij * sr12 v6(i,j) = scl(i,j)*epsij * sr6 vc(i,j) = scl(i,j)*qij / dij v12ij= v12ij + v12(i,j) v6ij = v6ij + v6(i,j) vcij = vcij + vc(i,j) else v12(i,j) = 0.00 v6(i,j) = 0.00 vc(i,j) = 0.00 endif totij= eb(i,j)+eth(i,j)+etr(i,j)+v12(i,j)-v6(i,j)+vc(i,j) write(iueij,'(x,i6,2(x,i3),7f9.3)') nconf,i,j,eb(i,j), : eth(i,j), etr(i,j), v12(i,j), -v6(i,j), vc(i,j), totij enddo enddo endif c calcula dipolo do soluto do i= 1, na(1) dipx = dipx + qi(i)*rx(i) dipy = dipy + qi(i)*ry(i) dipz = dipz + qi(i)*rz(i) enddo dip = sqrt(dipx*dipx + dipy*dipy + dipz*dipz) write(iudip,'(i6,4(2x,f8.4))') nconf, dipx*debye, dipy*debye, : dipz*debye, dip*debye te= teij + teijk + teijkl + v12ij - v6ij + vcij if (einit) write(iueint,'(x,i6,7f9.3)') nconf,teij,teijk, : teijkl,v12ij,-v6ij, vcij, te enddo 9999 continue close(iubnd) close(iuang) close(iudhd) close(iudip) close(iuxyz) close(iusol) if (finput .ne. 3) close(iuall) if (einit) then close(iuebnd) close(iueang) close(iuedhd) close(iueint) close(iueij) endif stop end c ******************************************************************* subroutine input(namein, namexyz, finput, ftinput, einit, ljname) implicit none common / atomic / n, na, nan, m, natot, amass, aicon integer mnt, mnat, mnats, tabe parameter (mnt= 6, mnat= 7000, mnats= 60, tabe= 105) integer n(mnt), na(mnt), nan(mnat), m, natot real amass(tabe) character*2 aicon(tabe) integer ieq, iflag, nini, nfim, finput, ftinput integer i, j, k, ii, iio, ki, kj character line*80 character namein*(*), namexyz*(*), ljname*(*), command*30 logical keyin, keyxyz, keymoltype, keynmol, keynatom, keylj logical flagstop, dtexst, einit c ******************************************************************* iflag = 0 finput= 1 ftinput= 1 einit = .false. keyin = .false. keyxyz = .false. keymoltype= .false. keynmol = .false. keynatom = .false. keylj = .false. flagstop = .false. do while (iflag .eq. 0) read(5,'(a)', end=999) line if ((index(line,'$end')) .ne. 0) then iflag= 1 else ieq = index(line,'nametop =') if(ieq .ne. 0) then keyin= .true. ieq= ieq + 9 read(line(ieq:),'(a)') command nini= 1 do while (command(nini:nini) .eq. ' ') nini= nini + 1 enddo nfim= nini do while (command(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo namein= command(nini:nfim) nfim= nfim - nini if ((nfim+1) .gt. 37) then write(6,*)'Error: File name is too large (Max.= 33', : 'character)' stop endif inquire(file=namein(:nfim), exist=dtexst) if (dtexst) then write(6,*) write(6,*)' Topology file of the solute: ##', : namein(:nfim),'##' else write(6,*)'Error: Topology file of the solute ##', : namein(:nfim), '## NOT found' stop endif write(6,*)' nametop = ', namein(:nfim) endif ieq = index(line,'nametraj =') if(ieq .ne. 0) then keyxyz= .true. ieq= ieq + 10 read(line(ieq:),'(a)') command nini= 1 do while (command(nini:nini) .eq. ' ') nini= nini + 1 enddo nfim= nini do while (command(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo namexyz= command(nini:nfim) nfim= nfim - nini if ((nfim+1) .gt. 37) then write(6,*)'Error: File name is too large (Max.= 33', : 'character)' stop endif inquire(file=namexyz(:nfim), exist=dtexst) if (dtexst) then write(6,*) write(6,*)' Trajectory file ##',namexyz(:nfim), '##' else write(6,*)'Error: Trajectory file ##', namexyz(:nfim), : '## NOT found' stop endif if ((index(namexyz,'.xyz')) .ne. 0) then ftinput= 2 write(6,*) write(6,*)' nametraj = ', namexyz(:nfim), ' (XYZ)' else write(6,*) write(6,*)' nametraj = ', namexyz(:nfim), ' (PDB)' endif endif ieq = index(line,'ljname =') if(ieq .ne. 0) then keylj= .true. ieq= ieq + 8 read(line(ieq:),'(a)') command nini= 1 do while (command(nini:nini) .eq. ' ') nini= nini + 1 enddo nfim= nini do while (command(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo ljname= command(nini:nfim) nfim= nfim - nini if ((nfim+1) .gt. 37) then write(6,*)'Error: File name is too large (Max.= 33', : 'character)' stop endif inquire(file=ljname(:nfim), exist=dtexst) if (dtexst) then write(6,*) write(6,*)' Dice ljname file ##',ljname(:nfim), '##' else write(6,*)'Error: Dice ljname file ##', ljname(:nfim), : '## NOT found' stop endif write(6,*)' LJname = ', ljname(:nfim) endif ieq = index(line,'format =') if(ieq .ne. 0) then ieq= ieq + 8 read(line(ieq:),*) finput if (finput .eq. 2) then write(6,*) write(6,*)' format = 2 (TINKER)' elseif (finput .eq. 1) then write(6,*) write(6,*)' format = 1 (GROMACS)' else write(6,*) write(6,*)' format = 3 (DICE)' finput= 3 endif endif endif ieq = index(line,'moltype =') if(ieq .ne. 0) then keymoltype= .true. ieq= ieq + 9 read(line(ieq:),*) m write(6,*)' moltype = ', m endif ieq = index(line,'nmol =') if(ieq .ne. 0) then keynmol= .true. ieq= ieq + 6 read(line(ieq:),*) (n(i),i=1,m) write(6,*)' nmol = ', (n(i),i=1,m) endif ieq = index(line,'natom =') if(ieq .ne. 0) then keynatom= .true. ieq= ieq + 7 read(line(ieq:),*) (na(i),i=1,m) write(6,*)' natom = ', (na(i),i=1,m) c read(5,*) (nan(i),i=1,na(1)) do i=1,na(1) read(5,*) nan(i) enddo iio= 1 ki = na(1) if (n(1).gt.1) then do j= 2, n(1) ki= ki+1 kj= ki+na(1)-1 ii= iio do i= ki, kj nan(i)= nan(ii) ii= ii + 1 enddo ki= kj enddo endif if (m .gt. 1) then do j= 2, m do k= 1, n(j) ki= ki+1 kj= ki+na(j)-1 if (k .eq. 1) then write(6,*) ki, kj read(5,*) (nan(i),i=ki,kj) write(6,*) (aicon(nan(i)),i=ki,kj) iio= ki ki = kj else ii= iio do i= ki, kj nan(i)= nan(ii) ii= ii + 1 enddo ki= kj endif enddo enddo endif natot = 0 do i= 1, m natot = natot + n(i)*na(i) enddo if (natot .gt. mnat) then write(6,*) 'Error: Increase "mnat" and recompiled' stop endif endif c ieq = index(line,'internal =') c if(ieq .ne. 0) then c if ((index(line,'yes').ne.0).or.(index(line,'YES').ne.0)) c : einit = .true. c endif enddo write(6,*) 'Total number of atoms: ', natot c if (einit) then c write(6,*) c write(6,*)' Internal energy will be calculate ' c else c write(6,*) c write(6,*)' Internal energy will NOT be calculate ' c endif write(6,*) write(6,*) if (.not. keyin) then write(6,*)'keyword "nametop" is missing in the input file' write(6,*)'nametop = name.key or name.top' flagstop= .true. endif if (.not. keyxyz) then write(6,*)'keyword "namexyz" is missing in the input file' write(6,*)'nametraj = name.xyz or name.pdb' flagstop= .true. endif if (.not. keymoltype) then write(6,*)'keyword "moltype" is missing in the input file' write(6,*)'moltype = number of differente molecules type' flagstop= .true. endif if (.not. keynmol) then write(6,*)'keyword "nmol" is missing in the input file' write(6,*)'nmol = quantity of molecules of each type' flagstop= .true. endif if (.not. keynatom) then write(6,*)'keyword "natom" is missing in the input file' write(6,*)'natom = quantity of atoms of each molecule type' flagstop= .true. endif if (finput .eq. 3) then if (.not. keylj) then write(6,*)'keyword "ljname" is missing in the input file' write(6,*)'ljname = molecule informations for Dice' flagstop= .true. endif endif if (flagstop) then stop endif 999 continue return end c ******************************************************************* C Funcao para calcular distancia entre os atomos i e j c ******************************************************************* function distancia(xi,yi,zi,xj,yj,zj) implicit none real distancia, xi, yi, zi, xj, yj, zj, xij, yij, zij xij = xi - xj yij = yi - yj zij = zi - zj distancia= sqrt(xij*xij + yij*yij + zij*zij) return end c ******************************************************************* C Funcao para calcular angulo entre os atomos i, j e k c ******************************************************************* function angulo(xi,yi,zi,xj,yj,zj,xk,yk,zk) implicit none c transforma radiano para graus real rad2deg parameter (rad2deg=57.29577951) real angulo, xi, yi, zi, xj, yj, zj, xk, yk, zk real xij, yij, zij, xkj, ykj, zkj, dij, dkj, escijk c distancia de j->i xij= xi-xj yij= yi-yj zij= zi-zj dij= sqrt(xij*xij + yij*yij + zij*zij) xij= xij / dij yij= yij / dij zij= zij / dij c distancia de j->k xkj= xk-xj ykj= yk-yj zkj= zk-zj dkj= sqrt(xkj*xkj + ykj*ykj + zkj*zkj) xkj= xkj / dkj ykj= ykj / dkj zkj= zkj / dkj c produto escalar de j->i com j->k normalizados escijk= (xij*xkj + yij*ykj + zij*zkj) angulo= acos(escijk)*rad2deg return end c ******************************************************************* C Funcao para calcular angulo de torcao entre os atomos i, j, k e l c ******************************************************************* function torcao(xi,yi,zi,xj,yj,zj,xk,yk,zk,xl,yl,zl) implicit none c transforma radiano para graus real rad2deg parameter (rad2deg=57.29577951) real torcao, xi, yi, zi, xj, yj, zj, xk, yk, zk, xl, yl, zl real xij, yij, zij, xkj, ykj, zkj, dij, dkj real xlk, ylk, zlk, xjk, yjk, zjk, dlk, djk real nxijk, nyijk, nzijk, nxjkl, nyjkl, nzjkl, escnn real nijk, njkl c distancia de j->i xij= xi - xj yij= yi - yj zij= zi - zj dij= sqrt(xij*xij + yij*yij + zij*zij) xij= xij / dij yij= yij / dij zij= zij / dij c distancia de j->k xkj= xk - xj ykj= yk - yj zkj= zk - zj dkj= sqrt(xkj*xkj + ykj*ykj + zkj*zkj) xkj= xkj / dkj ykj= ykj / dkj zkj= zkj / dkj c produto vetorial de j->k com j->i nxijk= (ykj*zij - zkj*yij) nyijk= (zkj*xij - xkj*zij) nzijk= (xkj*yij - ykj*xij) c normalizacao do nijk nijk = sqrt(nxijk*nxijk + nyijk*nyijk + nzijk*nzijk) nxijk = nxijk / nijk nyijk = nyijk / nijk nzijk = nzijk / nijk c distancia de k->l xlk= xl - xk ylk= yl - yk zlk= zl - zk dlk= sqrt(xlk*xlk + ylk*ylk + zlk*zlk) xlk= xlk / dlk ylk= ylk / dlk zlk= zlk / dlk c distancia de k->j xjk= -xkj yjk= -ykj zjk= -zkj djk= dkj c produto vetorial de k->l com k->j nxjkl= (ylk*zjk - zlk*yjk) nyjkl= (zlk*xjk - xlk*zjk) nzjkl= (xlk*yjk - ylk*xjk) c normalizacao do njkl njkl = sqrt(nxjkl*nxjkl + nyjkl*nyjkl + nzjkl*nzjkl) nxjkl= nxjkl / njkl nyjkl= nyjkl / njkl nzjkl= nzjkl / njkl c produto escalar entre os vetores normais normalizados:nijk e njkl escnn= (nxijk*nxjkl + nyijk*nyjkl + nzijk*nzjkl) torcao= acos(escnn)*rad2deg return end c ******************************************************************* C Funcao que verifica se o diedro esta orientado no sentido horario c ou anti-horario. c Desenvolvida por: Cedrick Miranda Mello c Implementada por: Evanildo Gomes Lacerda Jr. c ******************************************************************* logical function horario(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4) implicit none real rad2deg parameter (rad2deg=57.29577951) real x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4 real pa(3),pb(3),pc(3),pd(3) real panew(3), pdnew(3), pdnewr(3) real a(3,3), b(3,3), c(3,3), ba(3,3), cba(3,3) real ux,uy,uz,hcb,dcb, gama real da(3), dar(3) real dist_adnew, dist_adnewr real costheta, sintheta, cosphi, sinphi integer i,j c Aloca valores iniciais horario = .false. c Ponto A pa(1) = x1 pa(2) = y1 pa(3) = z1 c Ponto B pb(1) = x2 pb(2) = y2 pb(3) = z2 c Ponto C pc(1) = x3 pc(2) = y3 pc(3) = z3 c Ponto D pd(1) = x4 pd(2) = y4 pd(3) = z4 c Posiciona o segundo ponto B na origem do i= 1, 3 pa(i) = pa(i) - pb(i) pc(i) = pc(i) - pb(i) pd(i) = pd(i) - pb(i) pb(i) = pb(i) - pb(i) enddo c Definindo o vetor do eixo de rotacao BC ux = pc(1) uy = pc(2) uz = pc(3) dcb= sqrt(ux*ux+uy*uy+uz*uz) c Calcula angulo muito pequeno para rotacao arbitraria em torno de U gama = 0.01/rad2deg c Angulos de projecao com eixo x (theta) e z (phi) hcb = sqrt(ux*ux+uy*uy) costheta = ux/hcb sintheta = uy/hcb cosphi = uz/dcb sinphi = hcb/dcb c Matrizes de rotacao c A = rotaciona os eixos para projetar em x a(1,1) = costheta a(2,2) = costheta a(3,3) = 1 a(1,2) = sintheta a(2,1) = -sintheta a(1,3) = 0 a(2,3) = 0 a(3,1) = 0 a(3,2) = 0 c B = rotaciona os eixos para projetar em z b(1,1) = cosphi b(2,2) = 1 b(3,3) = cosphi b(1,3) = -sinphi b(3,1) = sinphi b(1,2) = 0 b(2,1) = 0 b(3,2) = 0 b(2,3) = 0 c C = rotaciona em z do angulo arbitrario c(1,1) = cos(gama) c(2,2) = cos(gama) c(3,3) = 1 c(1,2) = sin(gama) c(2,1) = -sin(gama) c(1,3) = 0 c(2,3) = 0 c(3,1) = 0 c(3,2) = 0 c Matriz BA de projecao dos eixos do i= 1, 3 do j= 1, 3 ba(i,j)= b(i,1)*a(1,j)+b(i,2)*a(2,j)+b(i,3)*a(3,j) enddo enddo c Matriz CBA completa do i= 1, 3 do j= 1, 3 cba(i,j)= c(i,1)*ba(1,j)+c(i,2)*ba(2,j)+c(i,3)*ba(3,j) enddo enddo c Novas coordenadas do i= 1, 3 panew(i)= ba(i,1)*pa(1)+ba(i,2)*pa(2)+ba(i,3)*pa(3) pdnew(i)= ba(i,1)*pd(1)+ba(i,2)*pd(2)+ba(i,3)*pd(3) pdnewr(i)= cba(i,1)*pd(1)+cba(i,2)*pd(2)+cba(i,3)*pd(3) enddo dist_adnew = 0 dist_adnewr= 0 do i= 1,3 da(i) = pdnew(i) - panew(i) dar(i)= pdnewr(i) - panew(i) dist_adnew = dist_adnew + da(i)*da(i) dist_adnewr= dist_adnewr + dar(i)*dar(i) enddo c Caso com a rotacao arbitraria no sentido horario aproxime o eixo c entao o sentido eh horario if (dist_adnew .lt. dist_adnewr) then horario = .true. endif return end block data atom implicit none common / atomic / n, na, nan, m, natot, amass, aicon integer mnt, mnat, mnats, tabe parameter (mnt= 6, mnat= 7000, mnats= 60, tabe= 105) integer n(mnt), na(mnt), nan(mnat), m, natot real amass(tabe) character*2 aicon(tabe) integer i data (amass(i), i= 1, tabe) : / 1.008, 4.003, 6.939, 9.012, 10.811, 12.011, 14.007, : 15.999, 18.998, 20.183, 22.989, 24.312, 26.982, 28.086, : 30.974, 32.064, 35.453, 39.948, 39.102, 40.080, 44.956, : 47.900, 50.942, 51.996, 54.938, 55.847, 58.933, 58.710, : 63.540, 65.370, 69.720, 72.590, 74.922, 78.960, 79.909, : 83.800, 85.470, 87.620, 88.905, 91.220, 92.906, 95.940, : 98.000, 101.070, 102.905, 106.400, 107.870, 112.400, 114.820, : 118.690, 121.750, 127.600, 126.904, 131.300, 132.905, 137.340, : 138.910, 140.120, 140.907, 144.240, 147.000, 150.350, 151.960, : 157.250, 158.924, 162.500, 164.930, 167.260, 168.934, 173.040, : 174.970, 178.490, 180.948, 183.850, 186.200, 190.200, 192.200, : 195.090, 196.967, 200.590, 204.370, 207.190, 208.980, 210.000, : 210.000, 222.000, 223.000, 226.000, 227.000, 232.038, 231.000, : 238.030, 237.000, 242.000, 243.000, 247.000, 247.000, 249.000, : 254.000, 253.000, 256.000, 254.000, 257.000, 0.000, 0.000 /, : (aicon(i), i=1, tabe) / ' H', 'He', 'Li', 'Be', ' B', :' C', ' N', ' O', ' F', 'Ne', 'Na', 'Mg', 'Al', 'Si', ' P', ' S', :'Cl', 'Ar', ' K', 'Ca', 'Sc', 'Ti', ' V', 'Cr', 'Mn', 'Fe', 'Co', :'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', :' Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', :'Sn', 'Sb', 'Te', ' I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', :'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', :'Hf', 'Ta', ' W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', :'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', ' U', 'Np', :'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lw', 'Xx', :'Bq'/ end