program configurations include 'order.prm' real ax, ay, az, rm, tmas, massa, lx, ly, lz, dens real xhist, yhist, yhist1, yhist2, yhist3 real rhmax, thmax, emax, rahb, mas(mnt) real delr, rcut, vol, const, first, inter, limit real mroo, mthe, mehb, mroh, nhbmed(3), nmed real*8 acm, acvol, one, achbm(3), achb(3) real*8 acsvhb(mnt) integer mathb, nfim, confprint, normcharge, inmp, bq integer ti, an, nn, nia, ii, i, j, k, nmp, nmc, ang, cm,fre integer nhist, iout, jj, ji, jf, iiflag, nch, ig98, ihb integer jmiao, jcmil, jdmil, jmil, jcent, jdec, junit integer npt, nptii, ijkl, titmp, mindex integer ntpo, nap, nd, ltopznd, lbotznd, ltopg98, lmidg98 integer l, ktmp, tk, tl, lmin, nas character out7*7,out6*6,out5*5,out4*4,out3*3,out2*2,out1*1 character mofile*30, atype*2 character*35 namein, ordfil, dstfil, nameou, xyzfil, angfil character*38 nordfil character*35 ppfil, ttfil, ccfil, llfil, molcasin, quantumtop character*35 fig98, fig98pc, fiall, hbfile, dimfil, hbfiles character*35 quantummid, quantumbot, grfile, enfile, fizndpc character*4 ord, dst, pt, xyz, angle character*5 hbtype character*1 amiao, acmil, admil, amil, acent, adec, aunit character*1 rule, fil(30), lin(80), aii, bii character fall*8, fg98*6, fg98pc*8, fzndpc*8, line*80 logical lpr, lprhb, lprhbs, dimer, lprdu, dtexst, gr, flex equivalence (namein,fil) equivalence (line,lin) c ******************************************************************* c nmp = numero de moleculas que serao impressas lidas no "molprint = nmp" c write(6,'('' '')') write(6,'(''*************************************************'')') write(6,'(''* Program ORDER '')') write(6,'(''* Last modification: 03/Set/2018 '')') write(6,'(''*************************************************'')') write(6,'('' '')') write(6,'('' '')') write(6,'('' Keywords: '')') write(6,'(''ljname = file_name_of_geometry_and_LJ_paramaters '')') write(6,'(''inname = file_name_of_XYZ_coordinates '')') write(6,'(''nmol = number_of_molecules (n(i),i=1,N_types) '')') write(6,'(''dens = density_value or 0 for_NPT_simulation '')') write(6,'(''cm = number_of_the_atom (spherical) '')') write(6,'('' 0 for the center-of-mass (spherical) '')') write(6,'('' -1 (minimum distance)'')') write(6,'(''freeze = atf1, atf2, atf3 '')') write(6,'('' where atf2 will be at the origin '')') write(6,'('' atf1 will be at the x axis '')') write(6,'('' atf3 will be at the xy plane '')') write(6,'(''flexible = yes or no '')') write(6,'(''irdf = yes or no '')') write(6,'(''molprint = number_molecules number_point_charges '')') write(6,'(''IF molprint = 1 n '')') write(6,'('' normcharge = number_that_divides_the_charges '')') write(6,'(''printconfig = yes (to print each file) or no '')') write(6,'(''printformat = 1 (MOLCAS) 2 (ZINDO) 3 (GAUSSIAN) '')') write(6,'(''topfile = filename_toplines_MOLCAS/ZINDO/GUASSIAN'')') write(6,'(''midfile = filename_midlines_GUASSIAN '')') write(6,'(''botfile = filename_bottomlines_ZINDO '')') write(6,'(''blank = type '')') write(6,'('' where type= number_of_molecules_type '')') write(6,'(''printdummy = yes or no '')') write(6,'(''printinterval = interval_configuration_to_analize'')') write(6,'(''angle = yes or no '')') write(6,'(''IF angle = yes '')') write(6,'('' atomsangle = at_1 at_2 at_3 (solute) '')') write(6,'('' at_s_1 at_s_2 at_s_3 (solvent) '')') write(6,'(''hbond = yes or no '')') write(6,'(''IF hbond = yes '')') write(6,'('' hbondcriteria = R_OO Theta_OOH Energy_ij '')') write(6,'('' soluteacceptor = number_of_atoms_that_accept_H'')') write(6,'('' solute_atoms_that_accept_H '')') write(6,'('' solventdonor = number_of_pairs_that_donnor_H'')') write(6,'('' solvent_H solvent_O '')') write(6,'('' solventacceptor =number_of_atoms_that_accept_H'')') write(6,'('' solvent_atoms_that_accept_H '')') write(6,'('' solutedonor = number_of_pairs_that_donnor_H'')') write(6,'('' solute_H solute_O '')') write(6,'(''hbondsolvent = yes '')') write(6,'(''IF hbondsolvent = yes '')') write(6,'('' shells = first inter limit '')') write(6,'('' '')') write(6,'('' '')') c variaveis default nn= 9999999 cm = 0 ang = 0 fre = 0 nmp = 0 inmp= 0 nmc = 0 dens= 0.00 lx= 0.00 ly= 0.00 lz= 0.00 bq = 0 lpr = .false. lprhb = .false. lprhbs= .false. dimer = .false. lprdu = .false. ord = '.or' dst = '.dst' pt = '.cf' xyz = '.xyz' angle= '.ang' fall = '_all.xyz' fg98 = '_g.gjf' fg98pc = '_gpc.gjf' fzndpc = '_zpc.dat' nhbsla = 0 nhbsld = 0 do j= 1, mnt acsvhb(j) = dble(0.0) svhb(j) = 0 nhbsva(j) = 0 nhbsvd(j) = 0 do i= 1, mhb ahba(i,j)= 0 ahbd(i,j)= 0 enddo enddo iout = 2 confprint= 0 delr = 0.1 acvol = dble(0.0) acm = dble(0.0) one = dble(1.0) do i= 1, 3 achb(i) = dble(0.0) achbm(i) = dble(0.0) enddo do i= 1, mna do j= 1, mnt nat(i,j)= 0 enddo enddo do i= 1,3 acroo(i)= dble(0.0) acthe(i)= dble(0.0) acehb(i)= dble(0.0) acroh(i)= dble(0.0) nachb(i)= dble(0.0) enddo c lendo input call input(cm, iout, nmp, nch, nmc, npt, dens, lpr, : lprdu, mofile, namein, quantumtop, normcharge, quantummid, : quantumbot, bq, gr, flex) inmp = nmp nptii= npt+7 nfim= 1 do while (mofile(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo nfim= nfim - 1 inquire(file=mofile, exist=dtexst) if (dtexst) then write(6,*)'Openning file:####',mofile(:nfim),'####' else write(6,*)'Error: File: ####', mofile(:nfim),'#### is missing' stop endif open(unit= umo, file= mofile(:nfim), status= 'old') read(umo,'(a)') rule if ( .not. ((rule .ne. '*') .or. (rule .ne. '+'))) then write(6,*)'Error: File ljname need definition of sigma rule' stop endif read(umo,*) m if (m .gt. mnt) then write(6,*) ' increase the parameter "mnt"' write(6,*) ' mnt= ', mnt, ' and your system is ', m stop endif do j= 1, m mas(j)= 0.00 dum(j)= 0 mtya(j)= 0 titmp = 1 read(umo,*) na(j) na(j)= na(j) + 1 if (na(j) .gt. mna) then write(6,*) ' increase the parameter "mna"' write(6,*) ' mna= ', mna, ' and your system is ', na(j) stop endif do i= 2, na(j) read(umo,*) ti, an, ax, ay, az,z(i,j),eps(i,j),sig(i,j) ano(i,j)= an if (an .eq. 104) dum(j)= dum(j) + 1 am(i,j)= amass(an) mas(j)= mas(j) + am(i,j) if (ti .le. titmp) then if (ti .eq. titmp) titmp= titmp+1 else write(6,*)'Error: In labels of RDF of molecule ', : 'type: ', j,' (number: ',titmp, ' is missing)' stop endif tya(i,j) = ti nat(ti,j)= nat(ti,j) + 1 if (ti .gt. mtya(j)) mtya(j)= ti enddo write(6,*) 'read ', na(j)-1, ' atoms of molecule type: ', j enddo close(unit= umo) c le do input: angulos e hbonds call input_2(dimer, ang, rhmax, thmax, emax, fre, first, inter, : limit,lprhb, lprhbs) if ((m.gt.2).and.(lprhb).and.(nmp.ne.0)) then write(6,*) 'Error: Configurations with a fixed number of', : ' Hbonds can the generate only for 1+N (solute+solvent)', : ' system' stop endif if ((ang.eq.1).and.(m .eq. 1)) then if (na(1) .eq. 3) then dimfil=namein(:npt) // '.dim' open(unit=udi,file=dimfil,status='unknown') write(udi,'('' Conf. mol1 mol2 Rcm rmin Theta1 Th :eta2 Phi'')') ppfil=namein(:npt) // 'p.xyz' open(unit= upp, file= ppfil, status='unknown') ttfil=namein(:npt) // 't.xyz' open(unit= utt, file= ttfil, status='unknown') ccfil=namein(:npt) // 'c.xyz' open(unit= ucc, file= ccfil, status='unknown') llfil=namein(:npt) // 'l.xyz' open(unit= ull, file= llfil, status='unknown') endif endif if (iout .eq. 1) then C MOLCAS molcasin = quantumtop elseif (iout .eq. 3) then C GAUSSIAN inquire(file=quantumtop, exist=dtexst) if (dtexst) then nfim= 1 do while (quantumtop(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo nfim= nfim - 1 write(6,*)'Openning file:####',quantumtop(:nfim),'####' open(unit= umo, file= quantumtop(:nfim), status= 'old') ltopg98=1 do while (nfim .ne. 0) read(umo,'(a)', end=8) topg98(ltopg98) ltopg98= ltopg98 + 1 enddo 8 continue ltopg98= ltopg98 - 1 close(umo) write(6,*) ltopg98, ' lines were read.' else ltopg98= 0 endif inquire(file=quantummid, exist=dtexst) if (dtexst) then nfim= 1 do while (quantummid(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo nfim= nfim - 1 write(6,*)'Openning file:####',quantummid(:nfim),'####' open(unit= umo, file= quantummid(:nfim), status= 'old') lmidg98=1 do while (nfim .ne. 0) read(umo,'(a)', end=9) midg98(lmidg98) lmidg98= lmidg98 + 1 enddo 9 continue lmidg98= lmidg98 - 1 write(6,*) lmidg98, ' lines were read.' close(umo) else lmidg98= 0 endif else C ZINDO inquire(file=quantumtop, exist=dtexst) if (dtexst) then nfim= 1 do while (quantumtop(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo nfim= nfim - 1 write(6,*)'Openning file:####',quantumtop(:nfim),'####' open(unit= umo, file= quantumtop(:nfim), status= 'old') ltopznd=1 do while (nfim .ne. 0) read(umo,'(a)',end=10) topznd(ltopznd) ltopznd= ltopznd + 1 enddo 10 continue ltopznd= ltopznd - 1 write(6,*) ltopznd, ' lines were read.' close(umo) else ltopznd= 0 write(6,*) ' None line was read.' endif inquire(file=quantumbot, exist=dtexst) if (dtexst) then nfim= 1 do while (quantumbot(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo nfim= nfim - 1 write(6,*)'Openning file:####',quantumbot(:nfim),'####' open(unit= umo, file= quantumbot(:nfim), status= 'old') lbotznd=1 do while (nfim .ne. 0) read(umo,'(a)',end=11) botznd(lbotznd) lbotznd= lbotznd + 1 enddo 11 continue lbotznd= lbotznd - 1 write(6,*) lbotznd, ' lines were read.' close(umo) else lbotznd= 0 write(6,*) ' None line was read.' endif endif c inicializacao das variaveis tmas= 0.0 nt = 0 k= 0 ndst= 0 nang= 0 nhist= 0 do 13 j= 1, 101 histdst(j)= 0 histcos(j)= 0 histang(j,1)= 0 histang(j,2)= 0 histang(j,3)= 0 13 continue do 20 j= 1, m tmas= tmas + mas(j)*n(j) nt = nt + n(j) do 19 i= 1, n(j) k= k + 1 ty(k)= j 19 continue do 20 i = 2, na(j) z(i,j) = z(i,j) * e eps(i,j)= 2.0 * sqrt(eps(i,j)) if (rule .eq. '*') then sig(i,j)= sqrt(sig(i,j)) else sig(i,j)= sig(i,j)/2.0 endif 20 continue massa= tmas*10.0/nav junit= 0 jdec = 0 jcent= 0 jmil = 0 jdmil= 0 jcmil= 0 jmiao= 0 hbfile = namein(:npt) // '.hbd' hbfiles= namein(:npt) // '_ss.hbd' grfile = namein(:npt) // '.gr' enfile = namein(:npt) // '.eij' fiall = namein(:npt) // fall fig98 = namein(:npt) // fg98 fig98pc= namein(:npt) // fg98pc fizndpc= namein(:npt) // fzndpc if ((nmp.eq.1).and.(nch.gt.0).and.(fre.eq.1)) then if (iout .eq. 3) then c GAUSSIAN write(6,'(/'' ASEC will be generated, see file '',a)')fig98pc elseif (iout .eq. 2) then c ZINDO write(6,'(/'' ASEC will be generated, see file '',a)')fizndpc endif elseif (fre.eq.0) then write(6,'(/'' ATTENTION: ASEC is generated only using the '', : ''keywork "freeze"''/)') endif c preparando para o calculo da RDFs ijkl = 0 if (gr) then write(6,'(/'' Energies that will be calculated, see file '', : a)') enfile write(6,'(/'' RDFs that will be calculated, see file '',a)') : grfile write(6,'('' 1 -> between CM of mol. 1 and CM of '', : ''others'')') write(6,'('' 2 -> between all atoms of mol. 1 '', : ''and nearest atom of others '')') do i= 1, m if (n(i) .eq. 1) then ii= i + 1 else ii= i endif if ( ii .le. m ) then do j= ii, m do k= 1, mtya(i) do l= 1, mtya(j) if (iindex(i,j,k,l) .eq. 0) then ijkl= ijkl + 1 iindex(i,j,k,l)= ijkl ktmp= 2 do while (k .ne. tya(ktmp,i)) ktmp= ktmp + 1 enddo tk= ano(ktmp,i) ktmp= 2 do while (l .ne. tya(ktmp,j)) ktmp= ktmp + 1 enddo tl= ano(ktmp,j) write(6,'(i5,'' -> between atom '',a,i5, : '' of mol. type '',i5,'' and '',a,i5, : '' of mol. type'',i5)') ijkl+2, aicon(tk), k, i, : aicon(tl), l, j if (i .eq. j) then iindex(i,j,l,k)= ijkl else iindex(j,i,l,k)= ijkl endif endif enddo enddo enddo endif enddo mindex= ijkl write(6,*) if (mindex .gt. nin) then write(6,*)'Error: Increase the parameter "nin" to ', mindex+1 stop endif else write(6,'(/'' RDFs that will NOT be calculated ''/)') endif if (iout .eq. 3) then c GAUSSIAN open(unit= ug98,file= fig98,status= 'unknown') if ((nmp.eq.1).and.(nch.gt.0).and.(fre.eq.1)) then open(unit= ugpc,file= fig98pc,status= 'unknown') if (ltopg98 .gt. 0) then do ig98=1, ltopg98 write(ugpc,'(a)') topg98(ig98) enddo endif endif elseif (iout .eq. 2) then c ZINDO if ((nmp.eq.1).and.(nch.gt.0).and.(fre.eq.1)) then open(unit= uzpc,file= fizndpc,status= 'unknown') if (ltopznd .gt. 0) then do ig98=1, ltopznd write(uzpc,'(a)') topznd(ig98) enddo endif endif endif open(unit= uall, file= fiall, status= 'unknown') dstfil=namein(:npt) // dst open(unit= udst, file= dstfil, status= 'unknown') if (dimer) then write(udst,'('' No Mol Rcm Energy '', : ''Rmin(1) Rmin(2) Theta Phi'')') elseif( ang .gt. 0 ) then write(udst,'('' No Mol Rcm Energy '', : ''Rmin(1) i j Rmin(2) i j Angle(nn) Angle(nv1)'', : ''Angle(nv2)'')') else write(udst,'('' No Mol Rcm Energy '', : '' Rmin(1) i j Rmin(2) i j'')') endif c Abre arquivo *.hbd if (lprhb) then open(unit= uhb, file= hbfile, status= 'unknown') write(uhb,'(''# Criteria '',3f7.3)') rhmax, thmax, emax write(uhb,'(''# No Hbond Mol R_OO Th_OOH'', :'' Energy Rcm R_OH Dip_A Dip_B Dip_AB '', :''Th_dip'')') endif c Abre arquivo *_ss.hbd if (lprhbs) then open(unit= uhbs, file= hbfiles, status= 'unknown') write(uhbs,'(''# Criteria '',3f7.3)') rhmax, thmax, emax write(uhbs,'(''# Shells '',2f7.3)') first, inter, limit write(uhbs,'(''# No Mol1 Hbond Rcm1 Mol2 Rcm2'', :'' R_OO Th_OOH Energy R_OH Type'')') endif nfim= 1 do while (namein(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo nfim= nfim - 1 inquire(file=namein, exist=dtexst) if (dtexst) then write(6,*)'Openning file:####',namein(:nfim),'####' else write(6,*)'Error: File: ####', namein(:nfim),'#### is missing' stop endif open(unit= uin, file= namein(:nfim), status= 'old') c Inicio do loop sobre as configuracoes do j= 1, nn c Zerando as HB do ihb= 1, mnm hbflag(ihb) = .false. enddo c definicao do nome do arquivo junit= junit + 1 if (junit .eq. 10) then junit= 0 jdec = jdec + 1 if (jdec .eq. 10) then jdec = 0 jcent= jcent + 1 if (jcent .eq. 10) then jcent = 0 jmil = jmil + 1 if (jmil .eq. 10) then jmil = 0 jdmil = jdmil + 1 if (jdmil .eq. 10) then jdmil = 0 jcmil = jcmil + 1 if (jcmil .eq. 10) then jcmil = 0 jmiao = jmiao + 1 endif endif endif endif endif endif c write(6,*) jmiao, jcmil, jdmil, jmil, jcent, jdec, junit amiao= char(jmiao+48) acmil= char(jcmil+48) admil= char(jdmil+48) amil = char(jmil+48) acent= char(jcent+48) adec = char(jdec+48) aunit= char(junit+48) if (nn .gt. 999999) then out7= amiao//acmil//admil//amil//acent//adec//aunit elseif (nn .gt. 99999) then out6= acmil // admil // amil // acent // adec // aunit out7= '0' // out6 elseif (nn .gt. 9999) then out5= admil // amil // acent // adec // aunit out7= '00' // out5 elseif (nn .gt. 999) then out4= amil // acent // adec // aunit out7= '000' // out4 elseif (nn .gt. 99) then out3= acent // adec // aunit out7= '0000' // out3 elseif ( nn .gt. 9) then out2= adec // aunit out7= '00000' // out2 else out1= aunit out7= '000000' // out1 endif nameou= namein(:npt) // out7 // pt ordfil= namein(:npt) // out7 // ord xyzfil= namein(:npt) // out7 // xyz c Leitura da configuracao nmp = inmp read(uin,*,end=999) nia read(uin,'(a)',end=999) line if ( mod(j, nmc) .eq. 0 ) then open(unit= uou, file=nameou, status= 'unknown') write(uou,*) nia write(uou,'(a)') line c leitura ou definicao do tamanho da caixa e Rcut do 25 i= 1, 50 if (lin(i) .eq. 'L') then jj= i + 3 read(line(jj:),*) lx, ly, lz endif 25 continue if (lx .lt. 1.0) then if (dens .gt. 0.001) then lx= ( massa / dens) ** (1.0/3.0) ly= lx lz= lx else stop ' Invalid value of density ' endif endif if (lx .lt. ly) then if (lx .lt. lz) then lmin= lx else lmin= lz endif else if (ly .lt. lz) then lmin= ly else lmin= lz endif endif rcut = lmin / 2.0 if ((delr*mbin).lt.rcut) then write(6,*)'Error: Increase the parameter "mbin" to ', : int(rcut/delr)+2 stop endif c leitura das coordenadas cartesianas do 40 i= 1, nia read(uin,'(a)',end=999) line jj= 1 do while (lin(jj) .eq. ' ') jj= jj + 1 enddo ji = jj do while (lin(jj) .ne. ' ') jj= jj + 1 enddo jf = jj - 1 atype = line(ji:jf) read(line(jj:),*) ax, ay, az write(uou,'(2x,a,3f15.5)') atype, ax, ay, az 40 continue close(unit= uou) confprint= confprint + 1 if((normcharge.eq.1).or. : ((normcharge.gt.1).and.(confprint.le.normcharge))) then write(6,'(''Analysing configuration '',i6,'' L = '', :3f9.4,'' Dens = '',f9.4)') j, lx, ly, lz, massa/(lx*ly*ly) open(unit= uou, file=nameou, status= 'unknown') nhist= nhist + 1 vol = lx * ly * lz acvol = acvol + dble(vol) acm = acm + one c le as configuracoes e reposiciona o CM call readconf(lx, ly, lz, mas, rm, cm) c calcula as RDFs e energias entre os diferentes tipos c de moleculas if (gr) then call rdf(enfile, rcut, delr, lx, ly, lz, rule, : sngl(acm),j, cm) else call energy(enfile, rcut, delr, lx, ly, lz, rule, : sngl(acm),j, cm) endif if ((ang .eq. 1) .and. (m .eq. 1)) : call neighbor(j, lx, ly, lz, dimer) c calcula lista de vizinhos call dist(rm) c analise das HB entre solvente-solvente if(lprhbs) then call hbondss(rhmax, thmax, emax, first, inter, limit, : j, lx, ly, lz, nch, lpr, out5, nhbmed) do jj= 1, 3 if (nhbmed(jj).gt.0.0) then achb(jj) = achb(jj) + one achbm(jj)= achbm(jj) + dble(nhbmed(jj)) endif enddo endif c analise das HB entre soluto-solvente if(lprhb) then call hbonds(rhmax, thmax, emax, cm, j, mathb, nmp) if (mathb .eq. 0) then close(unit= uou, status='delete') goto 888 endif if (inmp .eq. 0) nmp= mathb + 1 endif open(unit= uxyz, file= xyzfil, status= 'unknown') open(unit= uord, file= ordfil, status= 'unknown') if (iout .eq. 3) then C GAUSSIAN if (ltopg98 .gt. 0) then do ig98=1, ltopg98 write(ug98,'(a)') topg98(ig98) enddo endif elseif (iout .eq. 2) then c ZINDO if (ltopznd .gt. 0) then do ig98=1, ltopznd write(uord,'(a)') topznd(ig98) enddo endif endif c calcula quantas moleculas explicita que serao impressas c as cargas pontuais são incluidas na subrotina writeconf do i= 1, mnt svhb(i) = 0 enddo ntpo= n(1) i = 1 nd = tim(i) nap = 0 if (lprdu) then nap = nap + (na(nd)-1) else nap = nap + (na(nd)-1-dum(nd)) endif if (nmp .gt. 1) then i= i + 1 cKKKKKKKKKKK do while (i .lt. mnm) c se tiver HB if (lprhb) then if (hbflag(i)) then nd = tim(i) if (lprdu) then nas = na(nd)-1 else nas = na(nd)-1-dum(nd) endif nap = nap + nas svhb(nd)= svhb(nd) + 1 c write(6,*)'Achei HB tipo:',nd,i,nas, c : nap, nmp, mathb, svhb(nd), hbflag(i) if (svhb(nd) .eq. nmp) then i= mnm+1 else i= i + 1 endif else i= i + 1 endif else c se nao tiver HB k= hist(i) if (k .ne. 0) then nd = tim(k) if (lprdu) then nas = na(nd)-1 else nas = na(nd)-1-dum(nd) endif nap = nap + nas svhb(nd)= svhb(nd) + 1 c write(6,*)'Achei molecula tipo:',nd,i,nas, c : nap, nmp, svhb(nd) if (svhb(nd) .eq. nmp-1) then i= mnm+1 else i= i + 1 endif else i= i + 1 endif endif enddo endif c write(6,*)'Solute:',nap,nmp c Faz rotacoes para parar a molecula 1 if (fre .eq. 1) call project if (int(sngl(acm)) .eq. 1) then do i= 2, na(ty(1)) rx1(i)= rx(1,i) ry1(i)= ry(1,i) rz1(i)= rz(1,i) enddo endif call writeconf(nia, ii, nmp, nch, ang, nap, j, iout, : lprhb, mathb, molcasin, lprdu, lx, ly,lz, confprint, : normcharge, lmidg98, bq, fre, flex, lbotznd) do i = 2, m if (lprhb) then if (svhb(i).gt. 0) : write(6,*) 'Written explicitly ', nap, 'atoms. ', : ' One solute + ', svhb(i), : ' HB of solvent type ', i acsvhb(i) = acsvhb(i) + dble(svhb(i)) else if (svhb(i).gt. 0) : write(6,*) 'Written explicitly ', nap, 'atoms. ', : ' One solute + ', svhb(i), : ' molecules of solvent type ', i acsvhb(i) = acsvhb(i) + dble(svhb(i)) endif enddo if ((ii.ne.nt).or.(ii.ne.nmp)) then write(6,'(''Only '',i4,'' molecules were printed in file: : '',a)') ii, ordfil else write(6,'(i4,'' molecules were printed in file: : '',a)') ii, ordfil endif if (lpr) then if (lprhb) then if (ii .lt. 10) then aii= char(0+48) bii= char(ii+48) nordfil= ordfil(:nptii+3) //'-'// aii // bii elseif (ii .lt. 100) then aii= char(int(ii/10)+48) bii= char(mod(ii,10)+48) nordfil= ordfil(:nptii+3) // '-'// aii // bii else if (iout .eq. 3) then c GAUSSIAN nordfil= ordfil(:nptii) // '.com' elseif (iout .eq. 1) then c MOLCAS nordfil= ordfil(:nptii) // '.ord' else c ZINDO nordfil= ordfil(:nptii) // '.dat' endif endif else if (iout .eq. 3) then c GAUSSIAN nordfil= ordfil(:nptii) // '.com' elseif (iout .eq. 1) then c MOLCAS nordfil= ordfil(:nptii) // '.ord' else c ZINDO nordfil= ordfil(:nptii) // '.dat' endif endif open(unit= uordn, file= nordfil, status= 'unknown') rewind(uord) iiflag=0 if (iout .eq. 3) then c GAUSSIAN if (ltopg98 .gt. 0) then do ig98=1, ltopg98 write(uordn,'(a)') topg98(ig98) enddo endif endif do while (iiflag .eq. 0) read(uord,'(a)',end=39) line write(uordn,'(a)') line enddo 39 continue if (iout .eq. 3) then c GAUSSIAN write(uordn,*) write(uordn,*) elseif (iout .eq. 2) then c ZINDO if (lbotznd .gt. 0) then do ig98=1, lbotznd write(uordn,'(a)') botznd(ig98) enddo endif endif close(unit= uordn) close(unit= uord, status='delete') endif if (lpr) then close(unit= uxyz) if (iout .eq. 3) then close(unit= uord, status='delete') else close(unit= uord) endif else close(unit= uxyz, status='delete') close(unit= uord, status='delete') endif close(unit= uou, status='delete') else confprint= confprint - 1 goto 999 endif 888 continue else do 41 i= 1, nia read(uin,'(a)') line c read(uin,*) atype, ax, ay, az 41 continue endif enddo 999 continue if (iout .eq. 3) then c GAUSSIAN if ((nmp.eq.1).and.(nch.gt.0).and.(fre.eq.1)) then write(ugpc,*) write(ugpc,*) endif elseif (iout .eq. 2) then c ZINDO write(uzpc,'('' $END '')') write(uzpc,*) if (lbotznd .gt. 0) then do ig98=1, lbotznd write(uzpc,'(a)') botznd(ig98) enddo endif endif write(6,*) write(6,'(i6,'' configurations were analyzed'')') confprint write(6,*) if (confprint.lt.normcharge) then write(6,*) ' ########################################' write(6,*) ' ATENTION: The ASEC is WRONG ' write(6,*) ' The "normcharge" keyword is not correct.' write(6,*) ' It should be ', confprint write(6,*) ' ########################################' endif if (lprhb) then write(6,*) if (nhbsla .gt. 0) then do i= 1, nhbsla do j = 2, m write(6,'(''Over '',i6,'' configuration, there are '', : i6,'' Hbonds in the acceptor atom '',a,i3, : '' of the solvent type: '',i3,''. ('',f5.2, '' in average)'')') : nhist, ahba(i,j), aicon(ano(athbsl(i)+1,1)), athbsl(i), j, : (1.0*ahba(i,j))/(1.0*nhist) enddo enddo endif if (nhbsld .gt. 0) then do i= 1, nhbsld do j = 2, m write(6,'(''Over '',i6,'' configuration, there are '', : i6,'' Hbonds in the donor group ('', a, i3, '','', a, i3, : '') of the solvent type: '', i3,''. ('',f5.2, '' in average)'' : )') nhist, ahbd(i,j), aicon(ano(athbslo(i)+1,1)), athbslo(i), : aicon(ano(athbslh(i)+1,1)), athbslh(i), j, : (1.0*ahbd(i,j))/(1.0*nhist) enddo enddo endif else write(6,*) do j = 2, m if (sngl(acsvhb(j)).gt.0.01) : write(6,'(''Over '',i6,'' configuration, there are '', : f5.2,'' molecules in average of the solvent type: '',i3)') : nhist, sngl(acsvhb(j)/dble(1.0*nhist)), j enddo endif write(6,*) if ((ang .eq. 1) .and. (m .eq. 1)) then angfil = namein(:npt) // angle open(unit= uang, file= angfil, status= 'unknown') do 50 j= 1, 101 xhist= (j-0.5)*10.00 yhist= real(histdst(j)) / real(ndst) * 100.0 if (na(1) .eq. 3) then yhist1= real(histang(j,1)) / real(nang) * 100.0 else yhist1= real(histang(j,1)) / real(nang) * 100.0 yhist2= real(histang(j,2)) / real(nang) * 100.0 yhist3= real(histang(j,3)) / real(nang) * 100.0 endif if (xhist .lt. 181.0) then if (na(1) .eq. 3) then write(uang,*) xhist, yhist1 else write(uang,*) xhist, yhist1, yhist2, yhist3 endif endif 50 continue write(uang,*) write(uang,*)'# Cossine' do 70 j= 1, 101 xhist= (j-1)*0.02 yhist= real(histcos(j)) / real(nang) * 100.0 write(uang,*) xhist, yhist 70 continue close(unit= uang) endif c Fecha os arquivos de H bonds e escreve o sumario if (lprhb) close(unit= uhb) if (lprhbs) then write(6,*) write(6,*)' Average over the solvent-solvent H bonds: ' write(6,*)' ' do i= 1, 3 if (i.eq.3) then hbtype='BULK ' elseif (i .eq. 2) then hbtype='INTER' else hbtype='FIRST' endif nmed= sngl(achbm(i)/achb(i)) mroo= sngl(acroo(i)/nachb(i)) mthe= sngl(acthe(i)/nachb(i)) mehb= sngl(acehb(i)/nachb(i)) mroh= sngl(acroh(i)/nachb(i)) write(6,'(a5,x,f8.3,x,f8.4,x,f9.2,x,f8.4,x,f8.4)') : hbtype, nmed, mroo, mthe, mehb, mroh enddo close(unit= uhbs) endif c Normaliza as RDFs e escreve no arquivo grfile if (gr) then const = 2.0 / 3.0 * pi / sngl(acvol / acm) call distribution(grfile, delr,const,sngl(acm),rcut) endif close(unit= uall) close(unit= udst) if (dimer) then close(unit= ucc) close(unit= ull) close(unit= upp) close(unit= utt) endif c close(unit= uin) stop end block data atom include 'order.prm' 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 c ******************************************************************* subroutine dipole(j, dip1, dipj, dip1j, thdip) include 'order.prm' integer j, i, ti, k real dip1, dipj, dip1j, esc, thdip real dipx,dipy,dipz,dipx1,dipy1,dipz1,dipxs,dipys,dipzs c ******************************************************************* i= 1 ti= ty(i) dipx1= 0.00 dipy1= 0.00 dipz1= 0.00 do k= 2, na(ti) dipx1 = dipx1 + z(k,ti) * (rx(i,k) - rx(i,1)) dipy1 = dipy1 + z(k,ti) * (ry(i,k) - ry(i,1)) dipz1 = dipz1 + z(k,ti) * (rz(i,k) - rz(i,1)) enddo dip1 = sqrt(dipx1*dipx1 + dipy1*dipy1 + dipz1*dipz1)*debye i= j ti= ty(i) dipx= 0.00 dipy= 0.00 dipz= 0.00 do k= 2, na(ti) dipx = dipx + z(k,ti) * (rx(i,k) - rx(i,1)) dipy = dipy + z(k,ti) * (ry(i,k) - ry(i,1)) dipz = dipz + z(k,ti) * (rz(i,k) - rz(i,1)) enddo dipj= sqrt(dipx*dipx + dipy*dipy + dipz*dipz)*debye esc= (dipx1*dipx + dipy1*dipy + dipz1*dipz)*debye*debye thdip= acos(esc/dip1/dipj)*180.0/pi dipxs= dipx + dipx1 dipys= dipy + dipy1 dipzs= dipz + dipz1 dip1j= sqrt(dipxs*dipxs + dipys*dipys + dipzs*dipzs)*debye return end c ******************************************************************* c ** subroutine to normalize and write out the correlation ** c ** functions to unit ugr. ** c ******************************************************************* subroutine distribution(grfile, delr, const, iprop, rcut) include 'order.prm' integer i, j, k, l, ib, natom, an1, an2, ktmp, intj real nrml(mbin), delr, const, nideal, rlower, rupper, vo real rmid, coord, cnst, gr, iprop, rn, rni, rnj, rnatl, coo real rcut, rntj, grs(mntcm), coos(mntcm), coords(mntcm) character grfile*(*) c ******************************************************************* open(unit= ugr, file= grfile, status= 'unknown') c rn = real(nt) write(ugr,'('' # RDFs calculated with '',f10.1, : '' configurations.'')') iprop do i= 1, mbin rlower= real(i - 1) * delr rupper= rlower + delr nideal= 2.0 * const * (rupper**3 - rlower**3) nrml(i)= 1.00 / nideal / iprop enddo cnst= 6.0 * const vo = 2.0 / 3.0 * pi / const c write(ugr,'('' # RDF between atom cm 0 of mol. type 1 and '', :''cm 0 of mol. type all '')') coo = 0.00 coord= 0.00 do ib= 1, mbin rmid = (real(ib)-0.5) * delr if ((rmid .lt. rcut).or.((m .eq. 2).and.(n(2) .eq.1))) then gr= sngl(histcm(1,ib) * dble(nrml(ib) / rn / 2.0)) coo = coo + sngl(histcm(1,ib) * vo / rn / 2.0) coord= coord + cnst*rn*gr*rmid*rmid*delr c write(ugr,'(4f9.3)') rmid, gr, coord, coo write(ugr,'(3f12.3)') rmid, gr, coord endif enddo write(ugr,*) c write(ugr,'('' # RDF between atom all 0 of mol. type 1 and '', :''nearest 0 of mol. type all '')') if (m .le. 2) then coo = 0.00 coord= 0.00 do ib= 1, mbin rmid = (real(ib)-0.5) * delr if (rmid .lt. rcut) then gr= sngl(histcm(2,ib) * dble(nrml(ib) / rn / 2.0)) coo = coo + sngl(histcm(2,ib) * vo / rn / 2.0) coord= coord + cnst*rn*gr*rmid*rmid*delr c write(ugr,'(4f9.3)') rmid, gr, coord, coo write(ugr,'(3f12.3)') rmid, gr, coord endif enddo else coo = 0.00 coord= 0.00 do intj = 3, m+1 rntj = real(n(intj-1)) coos(intj) = 0.00 coords(intj)= 0.00 enddo do ib= 1, mbin rmid = (real(ib)-0.5) * delr if (rmid .lt. rcut) then gr= sngl(histcm(2,ib) * dble(nrml(ib) / rn / 2.0)) coo = coo + sngl(histcm(2,ib) * vo / rn / 2.0) coord= coord + cnst*rn*gr*rmid*rmid*delr do intj = 3, m+1 grs(intj)=sngl(histcm(intj,ib)*dble(nrml(ib)/rntj/2.0)) coos(intj)=coos(intj)+sngl(histcm(intj,ib)*vo/rntj/2.0) coords(intj)=coords(intj)+cnst*rntj*grs(intj)*rmid*rmid*delr enddo c write(ugr,'(4f9.3)') rmid, gr, coord, coo if (m.eq.3) then write(ugr,'(7f12.3)') rmid, gr, coord, : grs(3),coords(3),grs(4),coords(4) elseif (m.eq.4) then write(ugr,'(9f12.3)') rmid, gr, coord, : grs(3),coords(3),grs(4),coords(4),grs(5),coords(5) else write(ugr,*) rmid, gr, coord, : (grs(intj),coords(intj),intj=3,(m+1)) endif endif enddo endif write(ugr,*) c if (n(1) .gt. 1) then i= 1 j= 1 rni= real(n(i)) rnj= real(n(j)) do k= 1, mtya(i) do l= k, mtya(j) natom= nat(k,i) * nat(l,j) rnatl= nat(l,j) ktmp = 2 do while (k .ne. tya(ktmp,i)) ktmp = ktmp + 1 enddo an1 = ano(ktmp,i) ktmp = 2 do while (l .ne. tya(ktmp,j)) ktmp = ktmp + 1 enddo an2 = ano(ktmp,j) write(ugr,'('' # RDF between atom '',a,i5,'' of mol. type : '',i5,'' and '',a,i5,'' of mol. type'',i5)') :aicon(an1), k, i, aicon(an2), l, j coo = 0.00 coord= 0.00 do ib= 1, mbin rmid = (real(ib)-0.5) * delr if (rmid .lt. rcut) then if (k .ne. l) then gr= sngl(histij(iindex(i,j,k,l),ib) * dble(nrml(ib)/ :real(natom) / rni / rnj / 2.0)) coo = coo + sngl(histij(iindex(i,j,k,l),ib) * vo / :real(natom) / rni / rnj / 2.0) else gr= sngl(histij(iindex(i,j,k,l),ib) * dble(nrml(ib)/ :real(natom) / rni / rnj)) coo = coo + sngl(histij(iindex(i,j,k,l),ib) * vo / :real(natom) / rni / rnj) endif coord= coord + cnst*rnj*gr*rmid*rmid*delr c write(ugr,'(4f9.3)') rmid, gr, coord, coo write(ugr,'(3f12.3)') rmid, gr, coord endif enddo write(ugr,*) enddo enddo endif if (m .gt. 1) then do i= 1, m rni= real(n(i)) do j= 2, m rnj= real(n(j)) if ((i .eq. j) .and. (rni .gt. 1)) then do k= 1, mtya(i) do l= k, mtya(j) if (i .eq. 1) then natom= nat(k,i) * nat(l,j) rnatl= nat(l,j) ktmp = 2 do while (k .ne. tya(ktmp,i)) ktmp = ktmp + 1 enddo an1 = ano(ktmp,i) ktmp = 2 do while (l .ne. tya(ktmp,j)) ktmp = ktmp + 1 enddo an2 = ano(ktmp,j) write(ugr,'('' # RDF between atom '',a,i5,'' of m :ol. type '',i5,'' and '',a,i5,'' of mol. type'',i5)') :aicon(an1), k, i, aicon(an2), l, j else natom= nat(k,i) * nat(l,j) rnatl= nat(l,j) ktmp = 1 do while (k .ne. tya(ktmp,i)) ktmp = ktmp + 1 enddo an1 = ano(ktmp,i) ktmp = 1 do while (l .ne. tya(ktmp,j)) ktmp = ktmp + 1 enddo an2 = ano(ktmp,j) write(ugr,'('' # RDF between atom '',a,i5,'' of m :ol. type '',i5,'' and '',a,i5,'' of mol. type'',i5)') :aicon(an1), k, i, aicon(an2), l, j endif coo = 0.00 coord= 0.00 do ib= 1, mbin rmid = (real(ib)-0.5) * delr if (rmid .lt. rcut) then if (k .ne. l) then gr= sngl(histij(iindex(i,j,k,l),ib) * :dble(nrml(ib) / real(natom) / rni / rnj / 2.0)) coo = coo + sngl(histij(iindex(i,j,k,l),ib)* :vo/real(natom) / rni / rnj / 2.0) else gr= sngl(histij(iindex(i,j,k,l),ib) * :dble(nrml(ib) / real(natom) / rni / rnj)) coo = coo + (histij(iindex(i,j,k,l),ib) * vo / :real(natom) / rni / rnj) endif coord= coord + cnst*rnj*gr*rmid*rmid*delr c write(ugr,'(4f9.3)')rmid,gr,coord, coo write(ugr,'(3f12.3)')rmid,gr,coord endif enddo write(ugr,*) enddo enddo elseif (i .ne. j) then do k= 1, mtya(i) do l= 1, mtya(j) if (i .eq. 1) then natom= nat(k,i) * nat(l,j) rnatl= nat(l,j) ktmp = 2 do while (k .ne. tya(ktmp,i)) ktmp = ktmp + 1 enddo an1 = ano(ktmp,i) ktmp = 2 do while (l .ne. tya(ktmp,j)) ktmp = ktmp + 1 enddo an2 = ano(ktmp,j) write(ugr,'('' # RDF between atom '',a,i5,'' of m :ol. type '',i5,'' and '',a,i5,'' of mol. type'',i5)') :aicon(an1), k, i, aicon(an2), l, j else natom= nat(k,i) * nat(l,j) rnatl= nat(l,j) ktmp = 1 do while (k .ne. tya(ktmp,i)) ktmp = ktmp + 1 enddo an1 = ano(ktmp,i) ktmp = 1 do while (l .ne. tya(ktmp,j)) ktmp = ktmp + 1 enddo an2 = ano(ktmp,j) write(ugr,'('' # RDF between atom '',a,i5,'' of m :ol. type '',i5,'' and '',a,i5,'' of mol. type'',i5)') :aicon(an1), k, i, aicon(an2), l, j endif coo = 0.00 coord= 0.00 do ib= 1, mbin rmid= (real(ib)-0.5) * delr if (rmid .lt. rcut) then gr= sngl(histij(iindex(i,j,k,l),ib) * :dble(nrml(ib) / real(natom) / rni / rnj / 2.0)) coo = coo + sngl(histij(iindex(i,j,k,l),ib)*vo/ :real(natom) / rni / rnj / 2.0) coord= coord + cnst*rnj*gr*rmid*rmid*delr c write(ugr,'(4f9.3)')rmid,gr,coord, coo write(ugr,'(3f12.3)')rmid,gr,coord endif enddo write(ugr,*) enddo enddo endif enddo enddo endif close(unit= ugr) return end c ******************************************************************* c ** subroutine to accumulate the radial distribution functions and** c ** the pairwise energy. ** c ******************************************************************* subroutine rdf(enfile, rcut, delr, lx, ly, lz, rule, rconf, : conf, cm) include 'order.prm' integer i, j, k, l, ti, tj, ind, bkl, bij, indg, indl, tij, mtij integer tyai(mna),tyaj(mna), anoi(mna), anoj(mna), nv(men) integer conf, cm, tirm1, tjrm1, tirm2, tjrm2 real rxi(mna), ryi(mna), rzi(mna), lx, ly, lz, vol, rcut real rxj(mna), ryj(mna), rzj(mna), lmax, rm1, rm2 real rxkl, rykl, rzkl, rklsq, rkl, rijsq, dij, rconf real rxcm, rycm, rzcm, arxcm, arycm, arzcm, rminij, delr real v12ii, v6ii, vcii, v12kl, v6kl, vckl, ekl,skl, zkl real sr2, sr6, sr12, sklsq, rci, rc3i, rc9i, const real v12lrc, v12lrci, vlrcr, v12lrctij real v6lrc, v6lrci, vlrca, v6lrctij real rn, nl, eik, sik, sik2, sik6, sik12, rnterm real*8 two, zero, v12(men), v6(men), vc(men) character enfile*(*), rule*(*), tmp1*1, tmp2*1, tmpvar*450 c ******************************************************************* enrg(nt,nt)= 0.00 two= dble(2.0) zero= dble(0.0) vol = lx *ly *lz const = 2.0 / 3.0 * pi / vol if (lx .gt. ly) then lmax= lx else lmax= ly endif if (lmax .lt. lz) lmax= lz c open the output file of the energy k= 6 tmpvar(1:k) = 'NMOVE ' do i= 1, m tmp1= char(i+48) ti= ty(i) do j= i, m tmp2= char(j+48) tj= ty(j) indg= j indl= i-1 tij = indg-indl -(indl*(indl-1)/2) +indl*m tmpvar(k:k+9)=' U('//tmp1//'-'//tmp2//') ' k= k+10 v12(tij)= zero v6(tij) = zero vc(tij) = zero nv(tij) = 0 enddo enddo mtij= tij if (int(rconf) .eq. 1) then open(unit=uen, file=enfile, status='unknown') k=k-1 write(uen,*) tmpvar(:k) else open(unit=uen, file=enfile,access='append', status='unknown') endif c calcula a correcao de longo alcance c write(6,'(/''Long-Range Correction (LJ) added in the energy, c :'')') v12lrc= 0.0 v6lrc = 0.0 rci = 1.000 / rcut rc3i = rci * rci * rci rc9i = rc3i * rc3i * rc3i rn= real(nt) rnterm= (rn*rn - rn)/2.0 j= m do l= 1, m nl= real(n(l)) v12lrci= 0.0 v6lrci = 0.0 do i= 2, na(j) if (sig(i,j) .gt. 1.0e-05) then do k= 2, na(l) if (sig(k,l) .gt. 1.0e-05) then eik = eps(i,j) * eps(k,l) if (rule .eq. '*') then sik = sig(i,j) * sig(k,l) else sik = sig(i,j) + sig(k,l) endif sik2 = sik * sik sik6 = sik2 * sik2 * sik2 sik12= sik6 * sik6 vlrca = nl*rn * eik*sik12 / 3.0 vlrcr = nl*rn * eik*sik6 v12lrci= v12lrci + vlrca v6lrci = v6lrci + vlrcr endif enddo endif enddo v12lrci= const * v12lrci * rc9i v6lrci = const * v6lrci * rc3i c write(6,'(''between molecules type'',i3, c : '' and'', i3, '':'', f10.5, 2x,''(kcal/mol)'')') j, l, c : (v12lrci-v6lrci)/rn v12lrc = v12lrc + v12lrci v6lrc = v6lrc + v6lrci enddo v12lrctij= v12lrc / rnterm v6lrctij = v6lrc / rnterm c comeca o loop para calculo das RDFs e energias do i= 1, nt - 1 ti= ty(i) enrg(i,i)= 0.00 do k= 1, na(ti) rxi(k) = rx(i,k) ryi(k) = ry(i,k) rzi(k) = rz(i,k) tyai(k)= tya(k,ti) anoi(k)= ano(k,ti) enddo do j= i + 1, nt tj= ty(j) do l= 1, na(tj) rxj(l)= rx(j,l) ryj(l)= ry(j,l) rzj(l)= rz(j,l) tyaj(l)= tya(l,tj) anoj(l)= ano(l,tj) enddo rxcm = rxi(1) - rxj(1) rycm = ryi(1) - ryj(1) rzcm = rzi(1) - rzj(1) arxcm= anint(rxcm/lx) * lx arycm= anint(rycm/ly) * ly arzcm= anint(rzcm/lz) * lz rxcm = rxcm - arxcm rycm = rycm - arycm rzcm = rzcm - arzcm rijsq= rxcm*rxcm + rycm*rycm + rzcm*rzcm dij = sqrt(rijsq) if (i .eq. 1) then c acumula os hidtogramas das RDF_cm e MDDF bij = int(dij / delr) + 1 if (bij .lt. mbin+1) histcm(1,bij)= histcm(1,bij) + two if (cm .eq. 0) then c armazena na variavel DISTIJ a distancia entre CM de i e j distij(i,j)= dij distij(j,i)= dij elseif (cm .gt. 0) then c armazena na variavel DISTIJ a distancia entre o o atomo cm c da molecula I e o centro-de-massa das moleculas J rxcm= rx(i,cm+1) - rxj(1) rycm= ry(i,cm+1) - ryj(1) rzcm= rz(i,cm+1) - rzj(1) rxcm= rxcm - arxcm rycm= rycm - arycm rzcm= rzcm - arzcm dij= sqrt(rxcm*rxcm + rycm*rycm + rzcm*rzcm) distij(i,j)= dij distij(j,i)= dij endif endif rminij= 1000.00 v12kl= 0.00 v6kl = 0.00 vckl = 0.00 indg= tj indl= ti-1 tij = indg-indl -(indl*(indl-1)/2) +indl*m nv(tij)= nv(tij) + 1 rm1= 2.0*lmax*lmax rm2= 2.0*rm1 c comeca o loop sobre os atomos de i e j do k= 2, na(ti) if (tyai(k) .ne. 0) then do l= 2, na(tj) if (tyaj(l) .ne. 0) then ind= iindex(ti,tj,tyai(k),tyaj(l)) rxkl= rxi(k) - rxj(l) - arxcm rykl= ryi(k) - ryj(l) - arycm rzkl= rzi(k) - rzj(l) - arzcm rklsq= rxkl*rxkl + rykl*rykl + rzkl*rzkl rkl = sqrt(rklsq) if ((anoi(k) .ne. 104).and.(anoj(l) .ne. 104) .and. : (rkl .lt. 0.8)) then write(6,*)'Error: Overlap in a configuration ' write(6,*)' Distance: ', rkl write(6,*)' Between atom ',k-1,' of molecule ',i write(6,*)rxi(1), ryi(1), rzi(1) write(6,*)rxi(k)+rxi(1),ryi(k)+ryi(1),rzi(k)+rzi(1) write(6,*)' and atom ',l-1,' of molecule ',j write(6,*)rxj(1), ryj(1), rzi(1) write(6,*)rxj(l)+rxj(1),ryj(l)+ryj(1),rzj(l)+rzj(1) stop endif c identifica os dois atomos mais proximos entre 1 e j if (i .eq. 1) then if (rklsq .lt. rm1) then rm2= rm1 tirm2= tirm1 tjrm2= tjrm1 rm1= rklsq tirm1= ano(k,ti) tjrm1= ano(l,tj) elseif (rklsq .lt. rm2) then rm2= rklsq tirm2= ano(k,ti) tjrm2= ano(l,tj) endif endif if (rkl .lt. rminij) rminij= rkl c acumula o histograma da RDF bkl = int(rkl / delr) + 1 if (bkl.lt.mbin+1) histij(ind,bkl)=histij(ind,bkl)+two if (dij .lt. rcut) then ekl = eps(k,ti) * eps(l,tj) if (ekl .gt. 1.0e-05) then c calcula o termo Lennard-Jones da energia if (rule .eq. '*') then skl = sig(k,ti) * sig(l,tj) else skl = sig(k,ti) + sig(l,tj) endif if (skl .gt. 1.0e-05) then sklsq = skl * skl sr2 = sklsq / rklsq sr6 = sr2 * sr2 * sr2 sr12 = sr6 * sr6 v12ii= ekl * sr12 v6ii = ekl * sr6 v12kl= v12kl + v12ii v6kl = v6kl + v6ii endif endif c calcula o termo Coulombiano da energia zkl = z(k,ti) * z(l,tj) if (abs(zkl) .gt. 1.0e-05) then vcii = zkl / rkl vckl = vckl + vcii endif endif endif enddo endif enddo if (i .eq. 1) then bij = int(rminij / delr) + 1 if (bij .lt. mbin+1) then histcm(2,bij)= histcm(2,bij) + two if (m.gt.2) histcm(tj+1,bij)= histcm(tj+1,bij) + two endif c armazena o atomo de j mais proximo de 1 timin1(j)= tirm1 tjmin1(j)= tjrm1 rmin1(j) = sqrt(rm1) c armazena a menor distancia entre 1 e j if (cm .lt. 0) then distij(i,j)= rmin1(j) distij(j,i)= rmin1(j) endif c armazena o segundo atomo de j mais proximo de 1 timin2(j)= tirm2 tjmin2(j)= tjrm2 rmin2(j) = sqrt(rm2) c armazena a energia de 1 com j c enrg(1,j)= dble(v12kl) + v12lrctij - dble(v6kl) c : -v6lrctij + dble(vckl) endif c armazena a energia de i com j enrg(i,j)= dble(v12kl) + v12lrctij - dble(v6kl) : -v6lrctij + dble(vckl) enrg(j,i)= enrg(i,j) c armazena as energias entre todos os tipos de moleculas v12(tij)= v12(tij) + dble(v12kl) + v12lrctij v6(tij) = v6(tij) + dble(v6kl) + v6lrctij vc(tij) = vc(tij) + dble(vckl) enddo enddo c escreve as energias write(uen,*) conf,(sngl(v12(tij)-v6(tij)+vc(tij)),' ',tij=1,mtij) close(unit=uen) return end c ******************************************************************* c ** subroutine to calculate the pairwise energy. ** ** c ******************************************************************* subroutine energy(enfile, rcut, delr, lx, ly, lz, rule, rconf, : conf, cm) include 'order.prm' integer i, j, k, l, ti, tj, ind, indg, indl, tij, mtij c integer bkl, bij integer tyai(mna),tyaj(mna), anoi(mna), anoj(mna), nv(men) integer conf, cm, tirm1, tjrm1, tirm2, tjrm2 real rxi(mna), ryi(mna), rzi(mna), lx, ly, lz, vol, rcut real rxj(mna), ryj(mna), rzj(mna), lmax, rm1, rm2 real rxkl, rykl, rzkl, rklsq, rkl, rijsq, dij, rconf real rxcm, rycm, rzcm, arxcm, arycm, arzcm, rminij, delr real v12ii, v6ii, vcii, v12kl, v6kl, vckl, ekl,skl, zkl real sr2, sr6, sr12, sklsq, rci, rc3i, rc9i, const real v12lrc, v12lrci, vlrcr, v12lrctij real v6lrc, v6lrci, vlrca, v6lrctij real rn, nl, eik, sik, sik2, sik6, sik12, rnterm c real*8 two real*8 zero, v12(men), v6(men), vc(men) character enfile*(*), rule*(*), tmp1*1, tmp2*1, tmpvar*450 c ******************************************************************* enrg(nt,nt)= 0.00 c two= dble(2.0) c zero= dble(0.0) vol = lx *ly *lz const = 2.0 / 3.0 * pi / vol if (lx .gt. ly) then lmax= lx else lmax= ly endif if (lmax .lt. lz) lmax= lz c open the output file of the energy k= 6 tmpvar(1:k) = 'NMOVE ' do i= 1, m tmp1= char(i+48) ti= ty(i) do j= i, m tmp2= char(j+48) tj= ty(j) indg= j indl= i-1 tij = indg-indl -(indl*(indl-1)/2) +indl*m tmpvar(k:k+9)=' U('//tmp1//'-'//tmp2//') ' k= k+10 v12(tij)= zero v6(tij) = zero vc(tij) = zero nv(tij) = 0 enddo enddo mtij= tij if (int(rconf) .eq. 1) then open(unit=uen, file=enfile, status='unknown') k=k-1 write(uen,*) tmpvar(:k) else open(unit=uen, file=enfile,access='append', status='unknown') endif c calcula a correcao de longo alcance c write(6,'(/''Long-Range Correction (LJ) added in the energy, c :'')') v12lrc= 0.0 v6lrc = 0.0 rci = 1.000 / rcut rc3i = rci * rci * rci rc9i = rc3i * rc3i * rc3i rn= real(nt) rnterm= (rn*rn - rn)/2.0 j= m do l= 1, m nl= real(n(l)) v12lrci= 0.0 v6lrci = 0.0 do i= 2, na(j) if (sig(i,j) .gt. 1.0e-05) then do k= 2, na(l) if (sig(k,l) .gt. 1.0e-05) then eik = eps(i,j) * eps(k,l) if (rule .eq. '*') then sik = sig(i,j) * sig(k,l) else sik = sig(i,j) + sig(k,l) endif sik2 = sik * sik sik6 = sik2 * sik2 * sik2 sik12= sik6 * sik6 vlrca = nl*rn * eik*sik12 / 3.0 vlrcr = nl*rn * eik*sik6 v12lrci= v12lrci + vlrca v6lrci = v6lrci + vlrcr endif enddo endif enddo v12lrci= const * v12lrci * rc9i v6lrci = const * v6lrci * rc3i c write(6,'(''between molecules type'',i3, c : '' and'', i3, '':'', f10.5, 2x,''(kcal/mol)'')') j, l, c : (v12lrci-v6lrci)/rn v12lrc = v12lrc + v12lrci v6lrc = v6lrc + v6lrci enddo v12lrctij= v12lrc / rnterm v6lrctij = v6lrc / rnterm c comeca o loop para calculo das energias do i= 1, nt - 1 ti= ty(i) enrg(i,i)= 0.00 do k= 1, na(ti) rxi(k) = rx(i,k) ryi(k) = ry(i,k) rzi(k) = rz(i,k) tyai(k)= tya(k,ti) anoi(k)= ano(k,ti) enddo do j= i + 1, nt tj= ty(j) do l= 1, na(tj) rxj(l)= rx(j,l) ryj(l)= ry(j,l) rzj(l)= rz(j,l) tyaj(l)= tya(l,tj) anoj(l)= ano(l,tj) enddo rxcm = rxi(1) - rxj(1) rycm = ryi(1) - ryj(1) rzcm = rzi(1) - rzj(1) arxcm= anint(rxcm/lx) * lx arycm= anint(rycm/ly) * ly arzcm= anint(rzcm/lz) * lz rxcm = rxcm - arxcm rycm = rycm - arycm rzcm = rzcm - arzcm rijsq= rxcm*rxcm + rycm*rycm + rzcm*rzcm dij = sqrt(rijsq) if (i .eq. 1) then cc acumula os hidtogramas das RDF_cm e MDDF c bij = int(dij / delr) + 1 c if (bij .lt. mbin+1) histcm(1,bij)= histcm(1,bij) + two if (cm .eq. 0) then c armazena na variavel DISTIJ a distancia entre CM de i e j distij(i,j)= dij distij(j,i)= dij elseif (cm .gt. 0) then c armazena na variavel DISTIJ a distancia entre o atomo cm c da molecula I e o centro-de-massa das moleculas J rxcm= rx(i,cm+1) - rxj(1) rycm= ry(i,cm+1) - ryj(1) rzcm= rz(i,cm+1) - rzj(1) rxcm= rxcm - arxcm rycm= rycm - arycm rzcm= rzcm - arzcm dij= sqrt(rxcm*rxcm + rycm*rycm + rzcm*rzcm) distij(i,j)= dij distij(j,i)= dij endif endif rminij= 1000.00 v12kl= 0.00 v6kl = 0.00 vckl = 0.00 indg= tj indl= ti-1 tij = indg-indl -(indl*(indl-1)/2) +indl*m nv(tij)= nv(tij) + 1 rm1= 2.0*lmax*lmax rm2= 2.0*rm1 c comeca o loop sobre os atomos de i e j do k= 2, na(ti) if (tyai(k) .ne. 0) then do l= 2, na(tj) if (tyaj(l) .ne. 0) then ind= iindex(ti,tj,tyai(k),tyaj(l)) rxkl= rxi(k) - rxj(l) - arxcm rykl= ryi(k) - ryj(l) - arycm rzkl= rzi(k) - rzj(l) - arzcm rklsq= rxkl*rxkl + rykl*rykl + rzkl*rzkl rkl = sqrt(rklsq) if ((anoi(k) .ne. 104).and.(anoj(l) .ne. 104) .and. : (rkl .lt. 0.8)) then write(6,*)'Error: Overlap in a configuration ' write(6,*)' Distance: ', rkl write(6,*)' Between atom ',k-1,' of molecule ',i write(6,*)rxi(1), ryi(1), rzi(1) write(6,*)rxi(k)+rxi(1),ryi(k)+ryi(1),rzi(k)+rzi(1) write(6,*)' and atom ',l-1,' of molecule ',j write(6,*)rxj(1), ryj(1), rzi(1) write(6,*)rxj(l)+rxj(1),ryj(l)+ryj(1),rzj(l)+rzj(1) stop endif c identifica os dois atomos mais proximos entre 1 e j if (i .eq. 1) then if (rklsq .lt. rm1) then rm2= rm1 tirm2= tirm1 tjrm2= tjrm1 rm1= rklsq tirm1= ano(k,ti) tjrm1= ano(l,tj) elseif (rklsq .lt. rm2) then rm2= rklsq tirm2= ano(k,ti) tjrm2= ano(l,tj) endif endif if (rkl .lt. rminij) rminij= rkl cc acumula o histograma da RDF c bkl = int(rkl / delr) + 1 c if (bkl.lt.mbin+1) histij(ind,bkl)=histij(ind,bkl)+two if (dij .lt. rcut) then ekl = eps(k,ti) * eps(l,tj) if (ekl .gt. 1.0e-05) then c calcula o termo Lennard-Jones da energia if (rule .eq. '*') then skl = sig(k,ti) * sig(l,tj) else skl = sig(k,ti) + sig(l,tj) endif if (skl .gt. 1.0e-05) then sklsq = skl * skl sr2 = sklsq / rklsq sr6 = sr2 * sr2 * sr2 sr12 = sr6 * sr6 v12ii= ekl * sr12 v6ii = ekl * sr6 v12kl= v12kl + v12ii v6kl = v6kl + v6ii endif endif c calcula o termo Coulombiano da energia zkl = z(k,ti) * z(l,tj) if (abs(zkl) .gt. 1.0e-05) then vcii = zkl / rkl vckl = vckl + vcii endif endif endif enddo endif enddo if (i .eq. 1) then c bij = int(rminij / delr) + 1 c if (bij .lt. mbin+1) histcm(2,bij)= histcm(2,bij) + two c armazena o atomo de j mais proximo de 1 timin1(j)= tirm1 tjmin1(j)= tjrm1 rmin1(j) = sqrt(rm1) c armazena a menor distancia entre 1 e j if (cm .lt. 0) then distij(i,j)= rmin1(j) distij(j,i)= rmin1(j) endif c armazena o segundo atomo de j mais proximo de 1 timin2(j)= tirm2 tjmin2(j)= tjrm2 rmin2(j) = sqrt(rm2) c armazena a energia de 1 com j c enrg(1,j)= dble(v12kl) + v12lrctij - dble(v6kl) c : -v6lrctij + dble(vckl) endif c armazena a energia de i com j enrg(i,j)= dble(v12kl) + v12lrctij - dble(v6kl) : -v6lrctij + dble(vckl) enrg(j,i)= enrg(i,j) c armazena as energias entre todos os tipos de moleculas v12(tij)= v12(tij) + dble(v12kl) + v12lrctij v6(tij) = v6(tij) + dble(v6kl) + v6lrctij vc(tij) = vc(tij) + dble(vckl) enddo enddo c escreve as energias write(uen,*) conf,(sngl(v12(tij)-v6(tij)+vc(tij)),' ',tij=1,mtij) close(unit=uen) return end c ******************************************************************* c ** subroutine to identify H bonds between the solute and solvent.** ** c ******************************************************************* subroutine hbonds(rh, th, eh, cm, iconf, mathb, nmp) include 'order.prm' integer i, j, k, mathb, nmp, tj integer l1, l2, l3, cm, iconf real vx13, vy13, vz13, vx12, vy12, vz12, vx23, vy23, vz23 real v13, v12, v23, rh, th, eh, costh, theta real dip1, dipj, dip1j, thdip c ******************************************************************* mathb = 0 c analisa as ligacoes de HB do soluto como aceitador if (nhbsla .gt. 0) then do 100 i= 1, nhbsla l1= athbsl(i) + 1 do j= 2, nt tj = ty(j) do 90 k= 1, nhbsvd(tj) l2= athbsvh(k,tj) + 1 l3= athbsvo(k,tj) + 1 vx13 = rx(1,l1) - rx(j,l3) vy13 = ry(1,l1) - ry(j,l3) vz13 = rz(1,l1) - rz(j,l3) vx12 = rx(1,l1) - rx(j,l2) vy12 = ry(1,l1) - ry(j,l2) vz12 = rz(1,l1) - rz(j,l2) vx23 = rx(j,l2) - rx(j,l3) vy23 = ry(j,l2) - ry(j,l3) vz23 = rz(j,l2) - rz(j,l3) v13 = sqrt(vx13*vx13 + vy13*vy13 + vz13*vz13) v12 = sqrt(vx12*vx12 + vy12*vy12 + vz12*vz12) v23 = sqrt(vx23*vx23 + vy23*vy23 + vz23*vz23) costh= (vx23*vx13 + vy23*vy13 + vz23*vz13) / v23 / v13 theta= acos(costh)* 180.0/pi if (v13 .lt. rh) then c passou no criterio de distancia para formacao de HB c write(6,*) v13, theta, enrg(1,j) if ( (theta .lt. th) .and. (enrg(1,j) .lt. eh) ) then c passou no criterio de angulo e energia para c formacao de HB if (.not.(hbflag(j))) then c a molecula j passou pela primeira vez no criterio hbflag(j)= .true. if (nmp .eq. 0) then c TODAS as moleculas selecionadas pelo criterio de c HB serao impressas mathb= mathb + 1 ahba(i,tj) = ahba(i,tj) + 1 call dipole(j, dip1, dipj, dip1j, thdip) write(uhb,'(i8,''( '',i2,'' - '',i2,x,i2,'' )'', : i5,f8.4,f7.2,f9.4,5(x,f8.4),f8.2)') iconf, l1-1, : l3-1, l2-1, j, v13, theta, enrg(1,j),sqrt(ri(j)), : v12, dip1, dipj, dip1j, thdip if (mathb .gt. mhb) then write(6,*)' Error: increase the parameter "mhb"' write(6,*)' mhb= ',mhb, ' and your system is ', : mathb stop endif athb(mathb) = j elseif (mathb .lt. (nmp-1)) then c a molecula selecionada pelo criterio de HB sera c impressa mathb= mathb + 1 ahba(i,tj) = ahba(i,tj) + 1 call dipole(j, dip1, dipj, dip1j, thdip) write(uhb,'(i8,''( '',i2,'' - '',i2,x,i2,'' )'', : i5,f8.4,f7.2,f9.4,5(x,f8.4),f8.2)') iconf, l1-1, : l3-1, l2-1, j, v13, theta, enrg(1,j),sqrt(ri(j)), : v12, dip1, dipj, dip1j, thdip if (mathb .gt. mhb) then write(6,*)' Error: increase the parameter "mhb"' write(6,*)' mhb= ',mhb, ' and your system is ', : mathb stop endif athb(mathb) = j else c a molecula selecionada pelo criterio de HB NAO c sera impressa write(6,*)' Warning 1: The molecule ', j,' form', : ' HB but was not printed' write(6,'(i8,''( '',i2,'' - '',i2,x,i2,'' )'',i5, : f8.4,f7.2,f9.4,2(x,f8.4))') iconf, l1-1, l3-1, : l2-1, j, v13, theta, enrg(1,j), sqrt(ri(j)), v12 endif else c a molecula j ja passou no criterio anteriormente write(6,*)' Warning 2: The molecule ', j,' form ', : 'HB with more than one atom' write(6,'(i8,''( '',i2,'' - '',i2,x,i2,'' )'',i5, : f8.4,f7.2,f9.4,2(x,f8.4))') iconf, l1-1, l3-1, : l2-1, j, v13, theta, enrg(1,j), sqrt(ri(j)), v12 endif else c nao passou no criterio de angulo e energia para c formacao de HB, por isso escreve um WARNING. if ((theta .lt. 41.00) .and. (enrg(1,j) .lt. eh)) : write(6,*)' Warning 3: In configuration ', iconf, : ' The molecule ', j,'does not pass in the angular ', : 'and energetic criteria:', v13, theta, enrg(1,j) endif endif 90 continue enddo 100 continue endif c analisa as ligacoes de HB do soluto como doador if (nhbsld .gt. 0) then do 200 i= 1, nhbsld l2= athbslh(i) + 1 l3= athbslo(i) + 1 do j= 2, nt tj = ty(j) do 190 k= 1, nhbsva(tj) l1= athbsv(k,tj) + 1 vx13 = rx(j,l1) - rx(1,l3) vy13 = ry(j,l1) - ry(1,l3) vz13 = rz(j,l1) - rz(1,l3) vx12 = rx(j,l1) - rx(1,l2) vy12 = ry(j,l1) - ry(1,l2) vz12 = rz(j,l1) - rz(1,l2) vx23 = rx(1,l2) - rx(1,l3) vy23 = ry(1,l2) - ry(1,l3) vz23 = rz(1,l2) - rz(1,l3) v13 = sqrt(vx13*vx13 + vy13*vy13 + vz13*vz13) v12 = sqrt(vx12*vx12 + vy12*vy12 + vz12*vz12) v23 = sqrt(vx23*vx23 + vy23*vy23 + vz23*vz23) costh= (vx23*vx13 + vy23*vy13 + vz23*vz13) / v23 / v13 theta= acos(costh)* 180.0/pi if (v13 .lt. rh) then c passou no criterio de distancia para formacao de HB if ((theta .lt. th) .and. (enrg(1,j) .lt. eh)) then c passou no criterio de angulo e distancia para c formacao de HB if (.not.(hbflag(j))) then c a molecula j passou pela primeira vez no criterio hbflag(j)= .true. if (nmp .eq. 0) then c TODAS as moleculas selecionadas pelo criterio de c HB serao impressas mathb= mathb + 1 ahbd(i,tj) = ahbd(i,tj) + 1 call dipole(j, dip1, dipj, dip1j, thdip) write(uhb,'(i8,''( '',i2,'' - '',i2,x,i2,'' )'', : i5,f8.4,f7.2,f9.4,5(x,f8.4),f8.2)') iconf, l1-1, : l3-1, l2-1, j, v13, theta, enrg(1,j),sqrt(ri(j)), : v12, dip1, dipj, dip1j, thdip if (mathb .gt. mhb) then write(6,*)' Error: increase the parameter "mhb"' write(6,*)' mhb= ',mhb, ' and your system is ', : mathb stop endif athb(mathb) = j elseif (mathb .lt. (nmp-1)) then c a molecula selecionada pelo criterio de HB sera c impressa mathb= mathb + 1 ahbd(i,tj) = ahbd(i,tj) + 1 call dipole(j, dip1, dipj, dip1j, thdip) write(uhb,'(i8,''( '',i2,'' - '',i2,x,i2,'' )'', : i5,f8.4,f7.2,f9.4,5(x,f8.4),f8.2)') iconf, l1-1, : l3-1, l2-1, j, v13, theta, enrg(1,j),sqrt(ri(j)), : v12, dip1, dipj, dip1j, thdip if (mathb .gt. mhb) then write(6,*)'Error: increase the parameter "mhb"' write(6,*)' mhb= ',mhb, ' and your system is ', : mathb stop endif athb(mathb) = j else c a molecula selecionada pelo criterio de HB NAO c sera impressa write(6,*)' Warning 1: The molecule ', j,' form', : ' HB but was not printed' write(6,'(i8,''( '',i2,'' - '',i2,x,i2,'' )'',i5, : f8.4,f7.2,f9.4,2(x,f8.4))') iconf, l1-1, l3-1, : l2-1, j, v13, theta, enrg(1,j), sqrt(ri(j)), v12 endif else c a molecula j ja passou no criterio anteriormente write(6,*)' Warning 2: The molecule ', j,' form ', : 'HB with more than one atom' write(6,'(i8,''( '',i2,'' - '',i2,x,i2,'' )'',i5, : f8.4,f7.2,f9.4,2(x,f8.4))') iconf, l1-1, l3-1, : l2-1, j, v13, theta, enrg(1,j), sqrt(ri(j)), v12 endif else c nao passou no criterio de angulo e energia para c formacao de HB, por isso escreve um WARNING. if ((theta .lt. 41.00) .and. (enrg(1,j) .lt. eh)) : write(6,*)' Warning 3: In configuration ', iconf, : ' The molecule ', j,'does not pass in the angular ', : 'and energetic criteria:', v13, theta, enrg(1,j) endif endif 190 continue enddo 200 continue endif return end c ******************************************************************* c ** subroutine to identify H bonds between solvent molecules. ** ** c ******************************************************************* subroutine hbondss(rh, th, eh, first, inter, limit, iconf, lx, ly, : lz, nch, lpr, confii,nhbmed) include 'order.prm' integer i, j, k, l, ihbtp, ii, ti, tj, nathbs, li, lj, il integer l1, l2, l3, iconf, iihb, nch, inmol(3), inhb(3) real vx13, vy13, vz13, vx12, vy12, vz12, vx23, vy23, vz23 real v13, v12, v23, rh, th, eh, costh, theta, cmi, cmj real first, inter, limit, co(3,mna), lx, ly, lz, nhbmed(3) character*5 hbtype,molii, confii, hbty(3) character*1 aii, bii, cii, dii character hbt(3)*1, hbsfil*30, moltype*3 logical lpr c ******************************************************************* hbt(1)='f' hbt(2)='i' hbt(3)='b' hbty(1)= 'FIRST' hbty(2)= 'INTER' hbty(3)= 'BULK ' do i= 1, 3 inmol(i) = 0 inhb(i) = 0 nhbmed(i)= 0 enddo do i= 2, nt cmi= sqrt(ri(i)) if (cmi .lt. limit) then iihb=0 if (cmi.le.first) then inmol(1)= inmol(1) + 1 moltype= 'IN ' elseif (cmi.le.inter) then inmol(2)= inmol(2) + 1 moltype= 'INT' else inmol(3)= inmol(3) + 1 moltype= 'OUT' endif do j= 2, nt if (j .ne. i) then tj = ty(j) cmj= sqrt(ri(j)) c if (cmj .lt. limit) then c analisa as ligacoes de H entre i (como aceitador) e j if ((cmi.lt.first).and.(cmj.lt.first)) then hbtype= hbty(1) ihbtp= 1 elseif ((cmi.gt.inter).and.(cmj.gt.inter)) then hbtype= hbty(3) ihbtp= 3 else hbtype= hbty(2) ihbtp= 2 endif do l= 1, nhbsva(tj) l1= athbsv(l,tj) + 1 do k= 1, nhbsvd(tj) l2= athbsvh(k,tj) + 1 l3= athbsvo(k,tj) + 1 vx13 = rx(i,l1) - rx(j,l3) vy13 = ry(i,l1) - ry(j,l3) vz13 = rz(i,l1) - rz(j,l3) vx12 = rx(i,l1) - rx(j,l2) vy12 = ry(i,l1) - ry(j,l2) vz12 = rz(i,l1) - rz(j,l2) vx23 = rx(j,l2) - rx(j,l3) vy23 = ry(j,l2) - ry(j,l3) vz23 = rz(j,l2) - rz(j,l3) v13 = sqrt(vx13*vx13 + vy13*vy13 + vz13*vz13) v12 = sqrt(vx12*vx12 + vy12*vy12 + vz12*vz12) v23 = sqrt(vx23*vx23 + vy23*vy23 + vz23*vz23) costh= (vx23*vx13+vy23*vy13+vz23*vz13)/v23/v13 theta= acos(costh)* 180.0/pi if (v13 .lt. rh) then c passou no criterio de distancia para formacao de HB if ((theta .lt. th).and.(enrg(i,j) .lt. eh))then c passou no criterio de angulo e energia para c formacao de HB iihb= iihb + 1 if (cmi.le.first) then inhb(1)= inhb(1) + 1 elseif (cmi.le.inter) then inhb(2)= inhb(2) + 1 else inhb(3)= inhb(3) + 1 endif dii= char(iihb+48) write(uhbs,'(i8,x,i5,'' (A) '',x, f8.4,x,i5, : x,f8.4,f9.4,f7.2,2(x,f8.4),2x,a5,x,a3)') : iconf, i, cmi, j, cmj, v13, theta, enrg(i,j), : v12, hbtype, moltype acroo(ihbtp)= acroo(ihbtp) + dble(v13) acthe(ihbtp)= acthe(ihbtp) + dble(theta) acehb(ihbtp)= acehb(ihbtp) + dble(enrg(i,j)) acroh(ihbtp)= acroh(ihbtp) + dble(v12) nachb(ihbtp)= nachb(ihbtp) + dble(1.0) if (lpr) then if (i .lt. 10) then aii= char(0+48) bii= char(0+48) cii= char(i+48) molii= aii//bii//cii//'-'//dii elseif (i .lt. 100) then aii= char(0+48) bii= char(int(i/10)+48) cii= char(mod(i,10)+48) molii= aii//bii//cii//'-'//dii elseif (i .lt. 1000) then ii= int(i/100) aii= char(ii+48) ii= i - ii*100 bii= char(int(ii/10)+48) cii= char(mod(ii,10)+48) molii= aii//bii//cii//'-'//dii endif hbsfil= 'hbs'//hbt(ihbtp)//molii//'-'//confii//'.xyz' open(unit=usxyz,file=hbsfil,status='unknown') ti= tim(i) tj= tim(j) nathbs= na(ti) + na (tj) - 2 write(usxyz,*) nathbs write(usxyz,'('' Configuration number : '', : i7,''L = '',3f9.4)') iconf, lx, ly, lz do li= 2, na(ti) co(1,li)= rx(i,li) co(2,li)= ry(i,li) co(3,li)= rz(i,li) write(usxyz,'(2x,a,3f15.5)') : aicon(ano(li,ti)),(co(il,li),il=1,3) enddo do lj= 2, na(tj) co(1,lj)= rx(j,lj) co(2,lj)= ry(j,lj) co(3,lj)= rz(j,lj) write(usxyz,'(2x,a,3f15.5)') : aicon(ano(lj,tj)),(co(il,lj),il=1,3) enddo close(unit=usxyz) endif endif endif enddo enddo c analisa as ligacoes de H entre i (como doador) e j do l= 1, nhbsvd(tj) l2= athbsvh(l,tj) + 1 l3= athbsvo(l,tj) + 1 do k= 1, nhbsva(tj) l1= athbsv(k,tj) + 1 vx13 = rx(j,l1) - rx(i,l3) vy13 = ry(j,l1) - ry(i,l3) vz13 = rz(j,l1) - rz(i,l3) vx12 = rx(j,l1) - rx(i,l2) vy12 = ry(j,l1) - ry(i,l2) vz12 = rz(j,l1) - rz(i,l2) vx23 = rx(i,l2) - rx(i,l3) vy23 = ry(i,l2) - ry(i,l3) vz23 = rz(i,l2) - rz(i,l3) v13 = sqrt(vx13*vx13 + vy13*vy13 + vz13*vz13) v12 = sqrt(vx12*vx12 + vy12*vy12 + vz12*vz12) v23 = sqrt(vx23*vx23 + vy23*vy23 + vz23*vz23) costh= (vx23*vx13+vy23*vy13+vz23*vz13) /v23/v13 theta= acos(costh)* 180.0/pi if (v13 .lt. rh) then c passou no criterio de distancia para formacao de HB if ((theta .lt. th).and.(enrg(i,j) .lt. eh))then c passou no criterio de angulo e distancia para c formacao de HB iihb= iihb + 1 if (cmi.le.first) then inhb(1)= inhb(1) + 1 elseif (cmi.le.inter) then inhb(2)= inhb(2) + 1 else inhb(3)= inhb(3) + 1 endif dii= char(iihb+48) write(uhbs,'(i8,x,i5,'' (D) '',x, f8.4,x,i5, : x,f8.4,f9.4,f7.2,2(x,f8.4),2x,a5,x,a3)') : iconf, i, cmi, j, cmj, v13, theta, enrg(i,j), : v12, hbtype, moltype acroo(ihbtp)= acroo(ihbtp) + dble(v13) acthe(ihbtp)= acthe(ihbtp) + dble(theta) acehb(ihbtp)= acehb(ihbtp) + dble(enrg(i,j)) acroh(ihbtp)= acroh(ihbtp) + dble(v12) nachb(ihbtp)= nachb(ihbtp) + dble(1.0) if (lpr) then if (i .lt. 10) then aii= char(0+48) bii= char(0+48) cii= char(i+48) molii= aii//bii//cii//'-'//dii elseif (i .lt. 100) then aii= char(0+48) bii= char(int(i/10)+48) cii= char(mod(i,10)+48) molii= aii//bii//cii//'-'//dii elseif (i .lt. 1000) then ii= int(i/100) aii= char(ii+48) ii= i - ii*100 bii= char(int(ii/10)+48) cii= char(mod(ii,10)+48) molii= aii//bii//cii//'-'//dii endif hbsfil= 'hbs'//hbt(ihbtp)//molii//'-'//confii//'.xyz' open(unit=usxyz,file=hbsfil,status='unknown') ti= tim(i) tj= tim(j) nathbs= na(ti) + na (tj) - 2 write(usxyz,*) nathbs write(usxyz,'('' Configuration number : '', : i7,''L = '',3f9.4)') iconf, lx, ly, lz do li= 2, na(ti) co(1,li)= rx(i,li) co(2,li)= ry(i,li) co(3,li)= rz(i,li) write(usxyz,'(2x,a,3f15.5)') : aicon(ano(li,ti)),(co(il,li),il=1,3) enddo do lj= 2, na(tj) co(1,lj)= rx(j,lj) co(2,lj)= ry(j,lj) co(3,lj)= rz(j,lj) write(usxyz,'(2x,a,3f15.5)') : aicon(ano(lj,tj)),(co(il,lj),il=1,3) enddo close(unit=usxyz) endif endif endif enddo enddo c endif endif enddo endif enddo do i= 1, 3 if (inmol(i).ne. 0) nhbmed(i)= real(inhb(i))/real(inmol(i)) write(6,'(2a,i5,a,i5,a,f5.2)') hbty(i), ' Nmol= ', : inmol(i),' NHbond=', inhb(i), ' N_med= ', nhbmed(i) enddo return end c ******************************************************************* subroutine input(cm, iout, nmp, nch, nmc, npt, dens, lpr, : lprdu, mofile, namein, quantumtop, normcharge, quantummid, : quantumbot, bq, gr, flex) include 'order.prm' character mofile*(*),namein*(*),quantumtop*(*),quantummid*(*) character quantumbot*(*) character answp*4, answdu*4, answgr*4, answflex*4 logical lpr, lprdu, gr, flex integer cm, iout, nmp, nch, nmc, nfim, normcharge, bq real dens character line*80, lns(80)*1 integer i, is, ie, ieq, iflag, npt, nptii equivalence (line,lns) logical aljname, anamein, anmol, adens, acm, anmp, apconf logical aiout, atopfile, adummy, ainter, flagstop, abotfile logical amidfile, dtexst c ******************************************************************* iflag = 0 lpr = .false. lprdu = .false. gr = .false. flex = .false. aljname = .false. anamein = .false. anmol = .false. adens = .false. acm = .false. anmp = .false. apconf = .false. aiout = .false. atopfile= .false. amidfile= .false. abotfile= .false. adummy = .false. ainter = .false. flagstop= .false. normcharge = 1 open(unit= 7) do while (iflag .eq. 0) read(5,'(a)', end=999) line write(7,'(a)') line if ((index(line,'$end')) .ne. 0) then iflag= 1 else i= 1 ieq= 0 do while (i .lt. 80) if (lns(i) .eq. ' ') then i= i + 1 else is= i do while ((lns(i+1) .ne. ' ') .and. (lns(i+1) .ne. '=') : .and. ((i+1) .lt. 80)) i= i + 1 enddo if ((i+1) .lt. 80) then ie= i do while ((lns(i) .ne. '=') .and. (i .lt. 80)) i= i + 1 enddo if (i .lt. 80) then if (lns(i+1) .ne. ' ') then ieq= i else do while ((lns(i+1).eq.' ').and.((i+1).lt.80)) i= i + 1 enddo if ((i+1) .lt. 80) ieq= i + 1 endif i = 81 endif endif endif enddo if (ieq .ne. 0) then if((line(is:ie).eq.'ljname') .or. : (line(is:ie).eq.'LJNAME')) then read(line(ieq:),'(a)') mofile write(6,*)'read ljname= ',mofile aljname= .true. nfim= 1 do while (mofile(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo nfim= nfim - 1 inquire(file=mofile(:nfim), exist=dtexst) if (dtexst) then write(6,*)'Openning file:####',mofile(:nfim),'####' else write(6,*)'Error: File: ####', mofile(:nfim),'#### is' : ,' missing' stop endif open(unit= umo, file= mofile(:nfim), status= 'old') read(umo,*) read(umo,*) m close(umo) elseif((line(is:ie).eq.'inname') .or. : (line(is:ie).eq.'INNAME')) then read(line(ieq:),'(a)') namein write(6,*)'read inname= ', namein anamein= .true. nfim= 1 do while (namein(nfim:nfim) .ne. '.') nfim= nfim + 1 enddo npt= nfim - 1 nptii= npt+7 if ((nptii+1) .gt. 30) then write(6,*)'Error: File name is too large (Max.= 23 : character)' stop endif elseif((line(is:ie).eq.'nmol') .or. : (line(is:ie).eq.'NMOL')) then read(line(ieq:),*) (n(i),i=1,m) write(6,*)'read nmol= ', (n(i),i=1,m) anmol= .true. elseif((line(is:ie).eq.'dens') .or. : (line(is:ie).eq.'DENS'))then read(line(ieq:),*) dens write(6,*)'read dens= ', dens adens= .true. elseif((line(is:ie).eq.'cm') .or. : (line(is:ie).eq.'CM')) then read(line(ieq:),*) cm if (cm .gt. 0) then write(6,*)'read cm= ',cm,' (spherical with center in', : ' atom ', cm, ')' elseif (cm .lt. 0) then write(6,*)'read cm= ',cm,' (minimum distance) ' else write(6,*)'read cm= ',cm,' (spherical with center in', : ' the center-of-mass) ' endif acm= .true. elseif((line(is:ie).eq.'molprint') .or. : (line(is:ie).eq.'MOLPRINT')) then read(line(ieq:),*) nmp, nch write(6,*)'read molprint= ', nmp, nch anmp= .true. elseif((line(is:ie).eq.'printconfig') .or. : (line(is:ie).eq.'PRINTCONFIG')) then read(line(ieq:),'(a)') answp write(6,*)'read printconfig= ', answp apconf= .true. if ((answp(1:1).eq.'y').or.(answp(1:1).eq.'Y')) then lpr = .true. endif elseif((line(is:ie).eq.'printformat') .or. : (line(is:ie).eq.'PRINTFORMAT'))then read(line(ieq:),*) iout if (iout .eq. 1) then write(6,*)'read printformat= ',iout,' (MOLCAS)' elseif (iout .eq. 3) then write(6,*)'read printformat= ',iout,' (GAUSSIAN)' else write(6,*)'read printformat= ',iout,' (ZINDO)' endif aiout= .true. elseif((line(is:ie).eq.'blank') .or. : (line(is:ie).eq.'BLANK'))then read(line(ieq:),*) bq write(6,*)'read blank= ', bq elseif((line(is:ie).eq.'normcharge') .or. : (line(is:ie).eq.'NORMCHARGE'))then read(line(ieq:),*) normcharge write(6,*)'read normcharge= ', normcharge elseif((line(is:ie).eq.'topfile') .or. : (line(is:ie).eq.'TOPFILE')) then read(line(ieq:),'(a)') quantumtop write(6,*)'read topfile= ', quantumtop inquire(file=quantumtop, exist=dtexst) if (.not. dtexst) then write(6,*)' Error: topfile is missing' stop endif atopfile= .true. elseif((line(is:ie).eq.'midfile') .or. : (line(is:ie).eq.'MIDFILE')) then read(line(ieq:),'(a)') quantummid write(6,*)'read midfile= ', quantummid inquire(file=quantummid, exist=dtexst) if (.not. dtexst) then write(6,*)' Error: midfile is missing' stop endif abotfile= .true. elseif((line(is:ie).eq.'botfile') .or. : (line(is:ie).eq.'BOTFILE')) then read(line(ieq:),'(a)') quantumbot write(6,*)'read botfile= ', quantumbot inquire(file=quantumbot, exist=dtexst) if (.not. dtexst) then write(6,*)' Error: botfile is missing' stop endif abotfile= .true. elseif((line(is:ie).eq.'printdummy') .or. : (line(is:ie).eq.'PRINTDUMMY')) then read(line(ieq:),'(a)') answdu write(6,*)'read printdummy = ', answdu adummy= .true. if ((answdu(1:1).eq.'y').or.(answdu(1:1).eq.'Y')) : lprdu = .true. elseif((line(is:ie).eq.'printinterval') .or. : (line(is:ie).eq.'PRINTINTERVAL'))then read(line(ieq:),*) nmc write(6,*)'read printinterval= ', nmc ainter= .true. elseif((line(is:ie).eq.'irdf') .or. : (line(is:ie).eq.'IRDF')) then read(line(ieq:),'(a)') answgr write(6,*)'read irdf= ', answgr if ((answgr(1:1).eq.'y').or.(answgr(1:1).eq.'Y')) : gr = .true. elseif((line(is:ie).eq.'flexible') .or. : (line(is:ie).eq.'FLEXIBLE')) then read(line(ieq:),'(a)') answflex write(6,*)'read flexible= ', answflex if ((answflex(1:1).eq.'y').or.(answflex(1:1).eq.'Y')) : flex = .true. endif endif endif enddo 999 continue if (gr) then write(6,*)'irdf = yes (Will calculate the RDFcm-cm and MDDF)' else write(6,*)'irdf = no (Will not calculate the RDFcm-cm and ', : 'MDDF)' endif if (.not. aljname) then write(6,*)'keyword "ljname" is missing in the input file' write(6,*)'ljname = file_mane_of_geometry_and_LJparameters' flagstop= .true. endif if (.not. anamein) then write(6,*)'keyword "inname" is missing in the input file' write(6,*)'inname = file_name_of_XYZ_coordinates' flagstop= .true. endif if (.not. anmol) then write(6,*)'keyword "nmol" is missing in the input file' write(6,*)'nmol = number_of_molecules (n(i),i=1,N_types)' flagstop= .true. endif if (.not. adens) then write(6,*)'keyword "dens" is missing in the input file' write(6,*)'dens = density_value or 0 for_NPT_simulation' flagstop= .true. endif if (.not. acm) then write(6,*)'keyword "cm" is missing in the input file' write(6,*)'cm = number_of_the_atom (spherical)' write(6,*)' 0 for the center-of-mass (spherical)' write(6,*)' -1 (minimum distance)' flagstop= .true. endif if (.not. anmp) then write(6,*)'keyword "molprint" is missing in the input file' write(6,*)'molprint = number_molecules number_point_charges' flagstop= .true. endif if (.not. apconf) then write(6,*)'keyword "printconfig" is missing in the input file' write(6,*)'printconfig = yes (to print all files) or no ' flagstop= .true. endif if (.not. aiout) then write(6,*)'keyword "printformat" is missing in the input file' write(6,*)'printformat = 1 (MOLCAS) 2 (ZINDO) 3 (GAUSSIAN)' flagstop= .true. endif c if ((iout.eq.2).and.(.not. atopfile)) then c write(6,*)'keyword "topfile" is missing in the input file' c write(6,*)'topfile = file_name_toplines_ZINDO' c flagstop= .false. c endif c if ((iout.eq.2).and.(.not. abotfile)) then c write(6,*)'keyword "botfile" is missing in the input file' c write(6,*)'botfile = file_name_bottom_lines_ZINDO' c flagstop= .true. c endif if (.not. adummy) then write(6,*)'keyword "printdummy" is missing in the input file' write(6,*)'printdummy = yes or no ' flagstop= .true. endif if (.not. ainter) then write(6,*)'keyword "printinterval" is missing in the input', :' file' write(6,*)'printinterval = interval_configuration_to_analize' flagstop= .true. endif if (flagstop) then stop endif return end c ******************************************************************* subroutine input_2(dimer, ang, rhmax, thmax, emax, fre, first, : inter, limit, lprhb, lprhbs) include 'order.prm' character answ*4, answh*4, answhs*4 integer ang, fre real rhmax, thmax, emax, first, inter, limit character line*80, lns(80)*1 integer i, is, ie, ieq, iflag, j, at equivalence (line,lns) logical dimer, aangle, aatang, lprhb, lprhbs, ahbshell logical ahbond, ahbonds, acrit, asuac, asvdo, asudo, asvac c ******************************************************************* iflag = 0 ang = 0 fre = 0 lprhb = .false. lprhbs = .false. dimer = .false. aangle = .false. aatang = .false. ahbond = .false. ahbonds = .false. ahbshell= .false. acrit = .false. asuac = .false. asvdo = .false. asudo = .false. asvac = .false. rewind(7) do while (iflag .eq. 0) read(7,'(a)', end=999) line if ((index(line,'$end')) .ne. 0) then iflag= 1 else i= 1 ieq= 0 do while (i .lt. 80) if (lns(i) .eq. ' ') then i= i + 1 else is= i do while ((lns(i+1) .ne. ' ') .and. (lns(i+1) .ne. '=') : .and. ((i+1) .lt. 80)) i= i + 1 enddo if ((i+1) .lt. 80) then ie= i do while ((lns(i) .ne. '=') .and. (i .lt. 80)) i= i + 1 enddo if (i .lt. 80) then if (lns(i+1) .ne. ' ') then ieq= i else do while ((lns(i+1).eq.' ').and.((i+1).lt.80)) i= i + 1 enddo if ((i+1) .lt. 80) ieq= i + 1 endif i = 81 endif endif endif enddo if (ieq .ne. 0) then if((line(is:ie).eq.'angle') .or. : (line(is:ie).eq.'ANGLE'))then read(line(ieq:),'(a)') answ write(6,*)'read angle= ', answ aangle= .true. if ((answ(1:1).eq.'y').or.(answ(1:1).eq.'Y')) ang= 1 elseif((line(is:ie).eq.'atomsangle') .or. : (line(is:ie).eq.'ATOMSANGLE'))then read(line(ieq:),*) at1(1), at2(1), at3(1) write(6,*)'read atomsangle= ', at1(1), at2(1), at3(1) aatang= .true. if ((at1(1).eq.1).and.(at2(1).eq.1)) dimer=.true. if (m .gt. 1) then do j= 2, m read (7,*) at1(j), at2(j), at3(j) write(6,*) at1(j), at2(j), at3(j) enddo endif elseif((line(is:ie).eq.'freeze') .or. : (line(is:ie).eq.'FREEZE'))then read(line(ieq:),*) atf1, atf2, atf3 if ((atf2.gt.0).and.(atf2.lt.na(1))) then if ((atf1.gt.0).and.(atf1.lt.na(1))) then if ((atf3.gt.0).and.(atf3.lt.na(1))) then write(6,*)'read freeze= ', atf1, atf2, atf3 else if (na(1).eq.3) then write(6,*)'read freeze= ', atf1, atf2 write(6,*)'Used only for diatomic solute' write(6,*)'atf3 = 0' atf3 = 0 else write(6,*)'read freeze= ', atf1, atf2, atf3 write(6,*)'Error: atf3 = ', atf3, : ' is out of range' stop endif else if (na(1).eq.2) then write(6,*)'read freeze= ', atf1, atf2, atf3 write(6,*)'Used only for atomic solute' write(6,*)'atf1 = 0 and atf3 = 0' atf2 = 0 atf3 = 0 else write(6,*)'read freeze= ', atf1, atf2, atf3 write(6,*)'Error: atf1 = ', atf1, ' and atf3 = ', : atf3,' are out of range' stop endif endif c write(6,*)'read freeze= ', atf1, atf2, atf3 fre= 1 atf1= atf1 + 1 atf2= atf2 + 1 atf3= atf3 + 1 elseif((line(is:ie).eq.'hbond') .or. : (line(is:ie).eq.'HBOND'))then read(line(ieq:),'(a)') answh write(6,*)'read hbond= ', answh ahbond= .true. if ((answh(1:1).eq.'y').or.(answh(1:1).eq.'Y')) then lprhb= .true. if (n(1) .eq. 1) then else write(6,*)'Error: H bond is implemented only to :one solute and many solvent molecules' stop endif endif elseif((line(is:ie).eq.'hbondcriteria') .or. : (line(is:ie).eq.'HBONDCRITERIA'))then read(line(ieq:),*) rhmax, thmax, emax acrit= .true. write(6,*)'read hbondcriteria= ', rhmax, thmax, emax if (rhmax .lt. 1.0) rhmax= 4.001 if (thmax .lt. 1.0) thmax= 30.001 if (abs(emax) .lt. 0.5) emax= -0.001 elseif((line(is:ie).eq.'soluteacceptor') .or. : (line(is:ie).eq.'SOLUTEACCEPTOR'))then read(line(ieq:),*) nhbsla write(6,*)'read soluteacceptor= ', nhbsla asuac= .true. if (nhbsla .gt. mhb) then write(6,*) ' increase the parameter "mhb"' write(6,*) ' mhb= ', mhb, ' and your system is ', : nhbsla stop endif if (nhbsla .gt. 0 ) then read(7,*) (athbsl(i),i=1,nhbsla) do i=1,nhbsla write(6,'(x,a,i5)') aicon(ano(athbsl(i)+1,1)), : athbsl(i) enddo endif elseif((line(is:ie).eq.'solventdonor') .or. : (line(is:ie).eq.'SOLVENTDONOR'))then do j= 2, m if (j .gt. 2) read(7,'(a)', end=999) line write(6,*) ' Solvent type: ', j read(line(ieq:),*) nhbsvd(j) write(6,*)'read solventdonor= ', nhbsvd(j) asvdo= .true. if (nhbsvd(j) .gt. mhb) then write(6,*) ' Increase the parameter "mhb"' write(6,*) ' mhb= ', mhb, ' and your system is ', : nhbsvd(j) stop endif if (nhbsvd(j) .gt. 0 ) then do i= 1, nhbsvd(j) read(7,*) athbsvh(i,j), athbsvo(i,j) at= athbsvh(i,j)+1 if (ano(at,j) .ne. 1) then write(6,*)' Atom ', athbsvh(i,j), : ' is not H. It is a ', aicon(ano(at,j)) write(6,*)' ATTENTION: Write first the H number' endif write(6,*) aicon(ano(at,j)),athbsvh(i,j), : aicon(ano((athbsvo(i,j)+1),j)), athbsvo(i,j) enddo endif enddo elseif((line(is:ie).eq.'solutedonor') .or. : (line(is:ie).eq.'SOLUTEDONOR'))then read(line(ieq:),*) nhbsld write(6,*)'read solutedonor= ', nhbsld asudo= .true. if (nhbsld .gt. mhb) then write(6,*) ' increase the parameter "mhb"' write(6,*) ' mhb= ', mhb, ' and your system is ', : nhbsld stop endif if ((nhbsla+nhbsld) .eq. 0) then write(6,*)'Error: No HBond ?' stop endif if (nhbsld .gt. 0) then do i= 1, nhbsld read(7,*) athbslh(i), athbslo(i) at= athbslh(i)+1 if (ano(at,1) .ne. 1) then write(6,*)' Atom ', athbslh(i),' is not H. It is', : ' a ', aicon(ano(at,1)) write(6,*)' ATTENTION: Write first the H number' endif write(6,*) aicon(ano(at,1)), athbslh(i), : aicon(ano((athbslo(i)+1),1)), athbslo(i) enddo endif elseif((line(is:ie).eq.'solventacceptor') .or. : (line(is:ie).eq.'SOLVENTACCEPTOR'))then do j= 2, m if (j .gt. 2) read(7,'(a)', end=999) line write(6,*) ' Solvent type: ', j read(line(ieq:),*) nhbsva(j) write(6,*)'read solventacceptor= ', nhbsva(j) asvac= .true. if (nhbsva(j) .gt. mhb) then write(6,*) ' increase the parameter "mhb"' write(6,*) ' mhb= ', mhb, ' and your system is ', : nhbsva(j) stop endif if (nhbsva(j) .gt. 0) then read(7,*) (athbsv(i,j), i=1,nhbsva(j)) do i=1,nhbsva(j) write(6,'(x,a,i5)') aicon(ano(athbsv(i,j)+1,j)), : athbsv(i,j) enddo endif enddo elseif((line(is:ie).eq.'hbondsolvent') .or. : (line(is:ie).eq.'HBONDSOLVENT'))then read(line(ieq:),'(a)') answhs write(6,*)'read hbondsolvent= ', answhs ahbonds= .true. if ((answhs(1:1).eq.'y').or.(answhs(1:1).eq.'Y')) then lprhbs= .true. if (n(1) .eq. 1) then else write(6,*)'Error: H bond is implemented only to :one solute and many solvent molecules' stop endif endif elseif((line(is:ie).eq.'shells') .or. : (line(is:ie).eq.'SHELLS'))then read(line(ieq:),*) first, inter, limit write(6,*)'read shells= ', first, inter, limit ahbshell= .true. endif endif endif enddo 999 continue close(unit=7, status='delete') if (.not. aangle) then write(6,*)'keyword "angle" is missing in the input file' write(6,*)'angle = yes or no ' stop elseif ((ang.eq. 1) .and. (.not. aatang)) then write(6,*)'keyword "atomsangle" is missing in the input file' write(6,*)'atomsangle = at_1 at_2 at_3 (solute)' write(6,*)' at_s_1 at_s_2 at_s_3 (solvent)' stop endif if (.not. ahbond) then write(6,*)'keyword "hbond" is missing in the input file' write(6,*)'hbond = yes or no ' stop elseif (lprhb) then if (.not. acrit) then write(6,*)'keyword "hbondcriteria" is missing in the input', :' file' write(6,*)'hbondcriteria = R_OO Theta_OOH Energy_ij' stop endif if (.not. asuac) then write(6,*)'keyword "soluteacceptor" is missing in the input', :' file' write(6,*)'soluteacceptor = solute_atoms_that_accept_H' stop endif if (.not. asvdo) then write(6,*)'keyword "solventdonor" is missing in the input', :' file' write(6,*)'solventdonor = solvent_H solvent_O' stop endif if (.not. asvac) then write(6,*)'keyword "solventacceptor" is missing in the ', :'input file' write(6,*)'solventacceptor = solvent_atoms_that_accept_H' stop endif if (.not. asudo) then write(6,*)'keyword "solutedonor" is missing in the input', :' file' write(6,*)'solutedonor = solute_H solute_O' stop endif endif if ((ahbonds) .and. (.not. ahbond)) then write(6,*)'keyword "hbond" is missing in the input file' write(6,*)'hbond = yes or no ' stop elseif (lprhbs) then if (.not. ahbshell) then write(6,*)'keyword "shells" is missing in the input file' write(6,*)'shells = first inter limit' stop endif if (.not. acrit) then write(6,*)'keyword "hbondcriteria" is missing in the input', :' file' write(6,*)'hbondcriteria = R_OO Theta_OOH Energy_ij' stop endif if (.not. asuac) then write(6,*)'keyword "soluteacceptor" is missing in the input', :' file' write(6,*)'soluteacceptor = solute_atoms_that_accept_H' stop endif if (.not. asvdo) then write(6,*)'keyword "solventdonor" is missing in the input', :' file' write(6,*)'solventdonor = solvent_H solvent_O' stop endif if (.not. asvac) then write(6,*)'keyword "solventacceptor" is missing in the ', :'input file' write(6,*)'solventacceptor = solvent_atoms_that_accept_H' stop endif if (.not. asudo) then write(6,*)'keyword "solutedonor" is missing in the input', :' file' write(6,*)'solutedonor = solute_H solute_O' stop endif endif return end c ******************************************************************* subroutine neighbor(iconf, lx, ly, lz, dimer) include 'order.prm' integer i, j, k, irm, iconf, jo, ko integer l1, l2, l3 real ux12, uy12, uz12, ux13, uy13, uz13, unx, uny, unz, un real vx12, vy12, vz12, vx13, vy13, vz13, vnx, vny, vnz, vn real costh1, theta1, costh2, theta2, costh3, theta3 real x1, y1, z1, unc, v12, v13, d, lx, ly, lz, rxo, ryo, rzo real v1, v1x, v1y, v1z, v2, v2x, v2y, v2z real n1, n1x, n1y, n1z, n2, n2x, n2y, n2z logical dimer c ******************************************************************* c calcula a distribuicao angular para liquidos homogenios l1= at1(1) + 1 l2= at2(1) + 1 l3= at3(1) + 1 if (na(1) .eq. 3) then c Se liquidos diatomicos do 90 j= 1, n(1) k= irij(j) c vetor cm e a distancia entre os centros geometricos x1 = rx(k,1) - rx(j,1) y1 = ry(k,1) - ry(j,1) z1 = rz(k,1) - rz(j,1) x1= x1 - lx * anint(x1/lx) y1= y1 - ly * anint(y1/ly) z1= z1 - lz * anint(z1/lz) d= sqrt(x1*x1 + y1*y1 + z1*z1) if (dimer) then c para moleculas diatomicas homonucleores ja foram c identificados os dois atomos mais proximos entre as molec. c j e k. Aqui sao calculados os vetores v1, v2 e seu angulo. c v1= entre cm de j e o atomo mais distante. c v2= entre cm de k e o atomo mais distante. i= 5 - atjj(j) v1x = rx(j,i) - rx(j,1) v1y = ry(j,i) - ry(j,1) v1z = rz(j,i) - rz(j,1) v1 = sqrt(v1x*v1x + v1y*v1y +v1z*v1z) i= 5 - atij(j) v2x = rx(k,i) - rx(k,1) v2y = ry(k,i) - ry(k,1) v2z = rz(k,i) - rz(k,1) v2 = sqrt(v2x*v2x + v2y*v2y +v2z*v2z) theta1= acos((v1x*x1 + v1y*y1 + v1z*z1)/ v1 / d)* 180.0/pi theta2= acos((v2x*x1 + v2y*y1 + v2z*z1)/ v1 / d)* 180.0/pi if (theta1 .gt. 90.0) theta1= theta1 - 90.0 if (theta2 .gt. 90.0) theta2= theta2 - 90.0 c theta1= acos((v1x*v2x + v1y*v2y + v1z*v2z)/ v1 / v2)* 180.0/pi else c para moleculas diatomicas heteronucleores, os dois atomos c que definem os vetorees sao lidos inicialmente. Aqui sao c calculados os vetores v1, v2 e seu angulo. c v1= entre cm de j e o segundo atomo selecionados c v2= entre cm de k e o segundo atomo selecionados v1x = rx(j,l2) - rx(j,1) v1y = ry(j,l2) - ry(j,1) v1z = rz(j,l2) - rz(j,1) v1 = sqrt(v1x*v1x + v1y*v1y +v1z*v1z) v2x = rx(k,l2) - rx(k,1) v2y = ry(k,l2) - ry(k,1) v2z = rz(k,l2) - rz(k,1) v2 = sqrt(v2x*v2x + v2y*v2y +v2z*v2z) theta1= acos((v1x*x1 + v1y*y1 + v1z*z1)/ v1 / d)* 180.0/pi theta2= acos((v2x*x1 + v2y*y1 + v2z*z1)/ v1 / d)* 180.0/pi if (theta1 .gt. 90.0) theta1= theta1 - 90.0 if (theta2 .gt. 90.0) theta2= theta2 - 90.0 c theta1= acos((v1x*v2x + v1y*v2y + v1z*v2z)/ v1 / v2)* 180.0/pi endif c angulo dihedro n1x = v1y*z1 - v1z*y1 n1y = v1z*x1 - v1x*z1 n1z = v1x*y1 - v1y*x1 n1 = sqrt(n1x*n1x + n1y*n1y +n1z*n1z) n2x = v2y*z1 - v2z*y1 n2y = v2z*x1 - v2x*z1 n2z = v2x*y1 - v2y*x1 n2 = sqrt(n2x*n2x + n2y*n2y +n2z*n2z) theta3= acos((n1x*n2x + n1y*n2y + n1z*n2z)/ n1/n2)*180.0/pi if (theta3 .gt. 90.0) theta3= theta3 - 90.0 c acumulando no histograma nang= nang + 1 irm= int(theta1/10.0 + 1.0) histang(irm,1)= histang(irm,1) + 1 irm= int(theta2/10.0 + 1.0) histang(irm,2)= histang(irm,2) + 1 irm= int(theta3/10.0 + 1.0) histang(irm,3)= histang(irm,3) + 1 write(udi,'(i8,2i5,2f8.4,3f8.2)') iconf, j, k, d, : rij(j), theta1, theta2, theta3 jo = j ko = k if ( (theta1 .gt. 70.0) .and. (theta1 .lt. 90.0) .and. : (theta2 .gt. 70.0) .and. (theta2 .lt. 90.0) .and. : (theta3 .gt. 0.0) .and. (theta3 .lt. 20.0)) then write(upp,'('' 4 '')') write(upp,*) 'Configuracao', iconf,' Moleculas:', jo, ko i= 5 - atjj(j) rxo = rx(j,i) ryo = ry(j,i) rzo = rz(j,i) write(upp,'(a,3f12.6,2i3)')' O', rx(j,i)-rxo, ry(j,i)-ryo, : rz(j,i)-rzo i= atjj(j) write(upp,'(a,3f12.6,2i3)')' O', rx(j,i)-rxo, ry(j,i)-ryo, : rz(j,i)-rzo i= atij(j) k= irij(j) write(upp,'(a,3f12.6,2i3)')' O', rx(k,i)-rxo, ry(k,i)-ryo, : rz(k,i)-rzo i= 5 - atij(j) write(upp,'(a,3f12.6,2i3)')' O', rx(k,i)-rxo, ry(k,i)-ryo, : rz(k,i)-rzo endif if (((theta1 .gt.70.0) .and. (theta1 .lt. 90.0) .and. : (theta2 .gt. 0.0) .and. (theta2 .lt. 20.0)) .or. : ((theta1 .gt. 0.0) .and. (theta1 .lt. 20.0) .and. : (theta2 .gt.70.0) .and. (theta2 .lt. 90.0))) then write(utt,'('' 4 '')') write(utt,*)'Configuracao', iconf,' Moleculas:', jo, ko i= 5 - atjj(j) rxo = rx(j,i) ryo = ry(j,i) rzo = rz(j,i) write(utt,'(a,3f12.6,2i3)')' O', rx(j,i)-rxo, ry(j,i)-ryo, : rz(j,i)-rzo i= atjj(j) write(utt,'(a,3f12.6,2i3)')' O', rx(j,i)-rxo, ry(j,i)-ryo, : rz(j,i)-rzo i= atij(j) k= irij(j) write(utt,'(a,3f12.6,2i3)')' O', rx(k,i)-rxo, ry(k,i)-ryo, : rz(k,i)-rzo i= 5 - atij(j) write(utt,'(a,3f12.6,2i3)')' O', rx(k,i)-rxo, ry(k,i)-ryo, : rz(k,i)-rzo endif if ( (theta1 .gt.70.0) .and. (theta1 .lt. 90.0) .and. : (theta2 .gt.70.0) .and. (theta2 .lt. 90.0) .and. : (theta3 .gt.70.0) .and. (theta3 .lt. 90.0)) then write(ucc,'('' 4 '')') write(ucc,*)'Configuracao', iconf,' Moleculas:', jo, ko i= 5 - atjj(j) rxo = rx(j,i) ryo = ry(j,i) rzo = rz(j,i) write(ucc,'(a,3f12.6,2i3)')' O', rx(j,i)-rxo, ry(j,i)-ryo, : rz(j,i)-rzo i= atjj(j) write(ucc,'(a,3f12.6,2i3)')' O', rx(j,i)-rxo, ry(j,i)-ryo, : rz(j,i)-rzo i= atij(j) k= irij(j) write(ucc,'(a,3f12.6,2i3)')' O', rx(k,i)-rxo, ry(k,i)-ryo, : rz(k,i)-rzo i= 5 - atij(j) write(ucc,'(a,3f12.6,2i3)')' O', rx(k,i)-rxo, ry(k,i)-ryo, : rz(k,i)-rzo endif if ( (theta1 .gt. 0.0) .and. (theta1 .lt. 20.0) .and. : (theta2 .gt. 0.0) .and. (theta2 .lt. 20.0)) then write(ull,'('' 4 '')') write(ull,*)'Configuracao', iconf,' Moleculas:', jo, ko i= 5 - atjj(j) rxo = rx(j,i) ryo = ry(j,i) rzo = rz(j,i) write(ull,'(a,3f12.6,2i3)')' O', rx(j,i)-rxo, ry(j,i)-ryo, : rz(j,i)-rzo i= atjj(j) write(ull,'(a,3f12.6,2i3)')' O', rx(j,i)-rxo, ry(j,i)-ryo, : rz(j,i)-rzo i= atij(j) k= irij(j) write(ull,'(a,3f12.6,2i3)')' O', rx(k,i)-rxo, ry(k,i)-ryo, : rz(k,i)-rzo i= 5 - atij(j) write(ull,'(a,3f12.6,2i3)')' O', rx(k,i)-rxo, ry(k,i)-ryo, : rz(k,i)-rzo endif 90 continue else c Se liquidos poliatomicos do 100 j= 1, n(1) x1 = rx(j,l1) y1 = ry(j,l1) z1 = rz(j,l1) ux12 = rx(j,l2) - x1 uy12 = ry(j,l2) - y1 uz12 = rz(j,l2) - z1 ux13 = rx(j,l3) - x1 uy13 = ry(j,l3) - y1 uz13 = rz(j,l3) - z1 unx = uy12*uz13 - uy13*uz12 uny = uz12*ux13 - uz13*ux12 unz = ux12*uy13 - ux13*uy12 unc = unx*x1 + uny*y1 + unz*z1 un = sqrt(unx*unx + uny*uny + unz*unz) k= irij(j) x1 = rx(k,l1) y1 = ry(k,l1) z1 = rz(k,l1) d = unx*x1 + uny*y1 + unz*z1 if (d .lt. unc) then unx= -unx uny= -uny unz= -unz endif vx12 = rx(k,l2) - x1 vy12 = ry(k,l2) - y1 vz12 = rz(k,l2) - z1 vx13 = rx(k,l3) - x1 vy13 = ry(k,l3) - y1 vz13 = rz(k,l3) - z1 v12 = sqrt(vx12*vx12 + vy12*vy12 + vz12*vz12) v13 = sqrt(vx13*vx13 + vy13*vy13 + vz13*vz13) vnx = vy12*vz13 - vy13*vz12 vny = vz12*vx13 - vz13*vx12 vnz = vx12*vy13 - vx13*vy12 vn = sqrt(vnx*vnx + vny*vny + vnz*vnz) costh1= (unx*vnx + uny*vny + unz*vnz) / un / vn theta1= acos(costh1)* 180.0/pi costh2= (unx*vx12 + uny*vy12 + unz*vz12) / un / v12 theta2= acos(costh2)* 180.0/pi costh3= (unx*vx13 + uny*vy13 + unz*vz13) / un / v13 theta3= acos(costh3)* 180.0/pi c if (theta1 .gt. 89.999999) theta1= theta1 - 90.0 c if (theta2 .gt. 89.999999) theta2= theta2 - 90.0 c if (theta3 .gt. 89.999999) theta3= theta3 - 90.0 nang= nang + 1 irm= int((costh1+1.0)*50.0 + 1.0) c write(6,*) irm histcos(irm)= histcos(irm) + 1 irm= int(theta1/10.0 + 1.0) histang(irm,1)= histang(irm,1) + 1 irm= int(theta2/10.0 + 1.0) histang(irm,2)= histang(irm,2) + 1 irm= int(theta3/10.0 + 1.0) histang(irm,3)= histang(irm,3) + 1 100 continue endif return end c ******************************************************************* subroutine dist(rm) include 'order.prm' real rm integer i, j, k c ******************************************************************* do i= 0, mnm hist(i)= 0 enddo c O vetor DISTIJ armazenada as distancias entre as moleculas c 1 e J. Caso CM=-1, esta distancia eh a menor entre qualquer c atomo de 1 e qualquer atomo de J; caso CM=0, esta distancia eh c a distancia entre os centro-de-massa; e caso CM=numero_posistivo c esta eh a distancia entre o atomo CM da molecula 1 e o centro-de c -massa da molecula J. c O vetor HIST a lista de vizinhos da molecula 1 hist(0)= 0 distij(0,0) = 0.00 hist(1)= 1 distij(1,1) = 0.00 do i= 2, nt hist(i)= i j= i do while (distij(1,hist(j)) .lt. distij(1,hist(j-1))) k= hist(j-1) hist(j-1)= hist(j) hist(j)= k j= j - 1 enddo enddo return end c ******************************************************************* subroutine project include 'order.prm' integer i, j, tj real*8 xo, yo, zo, thx, thy, thz real*8 nvx, nvy, nvz, nv, cthx, sthx, cthy, sthy, sthz c ******************************************************************* xo= rx(1,atf2) yo= ry(1,atf2) zo= rz(1,atf2) do j= 1, nt tj= ty(j) do i= 2, na(tj) rx(j,i)= rx(j,i) - xo ry(j,i)= ry(j,i) - yo rz(j,i)= rz(j,i) - zo enddo enddo if (atf1 .gt. 1) then nvx = rx(1,atf1) nvy = ry(1,atf1) nvz = rz(1,atf1) nv = dsqrt(nvx*nvx + nvy*nvy + nvz*nvz) thz = dacos(nvz / nv) sthz= dsin(thz) thx = dacos(nvx / (nv*sthz)) write(6,*) 'Theta= ', thx*180.0/pi if (nvy .lt. 0.00) thx = -thx cthx = dcos(thx) sthx = dsin(thx) do j= 1, nt tj= ty(j) do i= 2, na(tj) xo = rx(j,i) yo = ry(j,i) rx(j,i)= xo*cthx + yo*sthx ry(j,i)= -xo*sthx + yo*cthx enddo enddo nvx = rx(1,atf1) nvy = ry(1,atf1) nvz = rz(1,atf1) thx = dacos(nvx / nv) c write(6,*) thx*180.0/pi if (nvz .lt. 0.00) thx = -thx sthx = dsin(thx) cthx = dcos(thx) do j= 1, nt tj= ty(j) do i= 2, na(tj) xo = rx(j,i) zo = rz(j,i) rx(j,i)= xo*cthx + zo*sthx rz(j,i)= -xo*sthx + zo*cthx enddo enddo endif if (atf3 .gt. 1) then nvx = rx(1,atf3) nvy = ry(1,atf3) nvz = rz(1,atf3) nv = dsqrt(nvx*nvx + nvy*nvy + nvz*nvz) thx = dacos(nvx / nv) sthx= dsin(thx) thy = dacos(nvy / (nv*sthx)) c write(6,*) thy*180.0/pi if (nvz .lt. 0.00) thy = -thy cthy = dcos(thy) sthy = dsin(thy) do j= 1, nt tj= ty(j) do i= 2, na(tj) yo = ry(j,i) zo = rz(j,i) ry(j,i)= yo*cthy + zo*sthy rz(j,i)= -yo*sthy + zo*cthy enddo enddo endif return end c ******************************************************************* subroutine readconf(lx, ly, lz, mas, rm, cm) include 'order.prm' integer ij, i, j, k, l, in, cm, irm real co(3,mna), mas(mnt),rxcm,rycm,rzcm real xcm, ycm, zcm, rxi, ryi, rzi, lx, ly,lz, rxl, ryl, rzl real anix, aniy, aniz, rm, rxo, ryo, rzo, rijsq, rmim character*2 at c ******************************************************************* c Inicia a leitura do arquivo de configuracoes read(uou,*) read(uou,*) if (cm .le. 0) then c se nenhum atomo foi selecionada, o centro de massa da primeira c molecula sera colocado no centro da caixa xcm = 0.0 ycm = 0.0 zcm = 0.0 do 70 l= 2, na(1) read(uou,'(2x,a,3f15.5)') at, (co(ij,l),ij=1,3) write(6,*) at, (co(ij,l),ij=1,3), lx,ly,lz rxl= co(1,l) ryl= co(2,l) rzl= co(3,l) xcm= xcm + am(l,1) * rxl ycm= ycm + am(l,1) * ryl zcm= zcm + am(l,1) * rzl 70 continue rxo= xcm / mas(1) ryo= ycm / mas(1) rzo= zcm / mas(1) else c se algum atomo foi selecionado, ele sera colocado no centro da c caixa do 80 l= 2, cm+1 read(uou,'(2x,a,3f15.5)') at, (co(ij,l),ij=1,3) rxl= co(1,l) ryl= co(2,l) rzl= co(3,l) 80 continue rxo= rxl ryo= ryl rzo= rzl endif write(6,*) '#CM = ', rxo, ryo, rzo c com o centro da caixa definido, o arquivo sera relido e sera c aplicado o metodo das imagem para reorganizar a caixa rewind(uou) read(uou,*) read(uou,*) in= 0 rm= 0.00 do 100 j= 1, m do 100 k= 1, n(j) in= in + 1 tim(in) = j xcm = 0.0 ycm = 0.0 zcm = 0.0 do 99 l= 2, na(j) read(uou,'(2x,a,3f15.5)') at, (co(ij,l),ij=1,3) rxl= co(1,l) - rxo ryl= co(2,l) - ryo rzl= co(3,l) - rzo xcm= xcm + am(l,j) * rxl ycm= ycm + am(l,j) * ryl zcm= zcm + am(l,j) * rzl rx(in,l)= rxl ry(in,l)= ryl rz(in,l)= rzl 99 continue rxi= xcm / mas(j) ryi= ycm / mas(j) rzi= zcm / mas(j) anix= anint(rxi/lx) aniy= anint(ryi/ly) aniz= anint(rzi/lz) rxi= rxi - lx * anix ryi= ryi - ly * aniy rzi= rzi - lz * aniz rx(in,1)= rxi ry(in,1)= ryi rz(in,1)= rzi ri(in)=rxi*rxi + ryi*ryi + rzi*rzi if (rm .lt. ri(in)) rm=ri(in) do 98 l= 2, na(j) rx(in,l)= rx(in,l) - lx * anix ry(in,l)= ry(in,l) - ly * aniy rz(in,l)= rz(in,l) - lz * aniz 98 continue 100 continue if ((m .gt. 1) .and. (n(1) .gt. 1)) then c se o soluto for composto de mais de uma molecula, o centro da c caixa sera definido como o centro de massa dentre todas as c moleculas do soluto j= 1 in= 0 rm= 0.00 xcm = 0.0 ycm = 0.0 zcm = 0.0 do 110 k= 1, n(j) xcm= xcm + rx(k,1) ycm= ycm + ry(k,1) zcm= zcm + rz(k,1) 110 continue rxo= xcm / real(n(j)) ryo= ycm / real(n(j)) rzo= zcm / real(n(j)) do 120 j= 1, m do 120 k= 1, n(j) in= in + 1 tim(in) = j xcm = 0.0 ycm = 0.0 zcm = 0.0 do 119 l= 2, na(j) rxl= rx(in,l) - rxo ryl= ry(in,l) - ryo rzl= rz(in,l) - rzo xcm= xcm + am(l,j) * rxl ycm= ycm + am(l,j) * ryl zcm= zcm + am(l,j) * rzl rx(in,l)= rxl ry(in,l)= ryl rz(in,l)= rzl 119 continue rxi= xcm / mas(j) ryi= ycm / mas(j) rzi= zcm / mas(j) anix= anint(rxi/lx) aniy= anint(ryi/ly) aniz= anint(rzi/lz) rxi= rxi - lx * anix ryi= ryi - ly * aniy rzi= rzi - lz * aniz rx(in,1)= rxi ry(in,1)= ryi rz(in,1)= rzi ri(in)=rxi*rxi + ryi*ryi + rzi*rzi if (rm .lt. ri(in)) rm=ri(in) do 118 l= 2, na(j) rx(in,l)= rx(in,l) - lx * anix ry(in,l)= ry(in,l) - ly * aniy rz(in,l)= rz(in,l) - lz * aniz 118 continue 120 continue elseif (m .eq. 1) then c caso o liquido sera homogenio, ou seja, composto por apenas c um tipo de molecula, sera identificada a molecula vizinha c mais proxima no vetor IRIJ. Adicionalmente, sera computado o c histograma destas distancias no vetor HISTDST do j= 1, n(1) irij(j)= 0 atij(j)= 0 atjj(j)= 0 rij(j) = 0.00 enddo do j= 1, n(1) rmim= 999.99 do k= 1, n(1) if (j .ne. k) then rxcm= rx(k,1) - rx(j,1) rycm= ry(k,1) - ry(j,1) rzcm= rz(k,1) - rz(j,1) anix= anint(rxcm/lx) aniy= anint(rycm/ly) aniz= anint(rzcm/lz) rxcm= rxcm - lx * anix rycm= rycm - ly * aniy rzcm= rzcm - lz * aniz if (cm .gt. 0) then c RIJ= distancia entre o atomo selecionado e CM rxi= rx(k,1) - rx(j,cm+1) ryi= ry(k,1) - ry(j,cm+1) rzi= rz(k,1) - rz(j,cm+1) rxi= rxi - lx * anix ryi= ryi - ly * aniy rzi= rzi - lz * aniz rijsq= rxi*rxi + ryi*ryi + rzi*rzi if (rmim .gt. rijsq) then rmim= rijsq rij(j) = sqrt(rijsq) irij(j)= k atij(j)= 1 atjj(j)= cm endif elseif (cm .lt. 0) then c RIJ= menor distancia entre atomos do i= 2, na(1) do l= 2, na(1) rxi= rx(k,l) - rx(j,i) ryi= ry(k,l) - ry(j,i) rzi= rz(k,l) - rz(j,i) rxi= rxi - lx * anix ryi= ryi - ly * aniy rzi= rzi - lz * aniz rijsq= rxi*rxi + ryi*ryi + rzi*rzi if (rmim .gt. rijsq) then rmim= rijsq rij(j) = sqrt(rijsq) irij(j)= k atij(j)= l atjj(j)= i endif enddo enddo else c RIJ= distancia entre os centros de massa rijsq= rxcm*rxcm + rycm*rycm + rzcm*rzcm if (rmim .gt. rijsq) then rmim= rijsq rij(j) = sqrt(rijsq) irij(j)= k atij(j)= 1 atjj(j)= 1 endif endif endif enddo irm= int(sqrt(rmim)*10.0 + 0.5) histdst(irm)= histdst(irm) + 1 ndst= ndst + 1 enddo endif return end c ******************************************************************* subroutine writeconf(nia, ii, nmp, nch, ang, nap, iconf, iout, : lprhb, mathb, molcasin, lprdu, lx, ly, lz, confprint, : normcharge, lmidg98, bq, fre, flex, lbotznd) include 'order.prm' integer nia, i, j, k, l, il, ii, mathb, nch, ta, nata, i0, jj integer nmp, nap, ang, l1, l2, l3, iconf, iout, nfim, nd, bq integer confprint, normcharge, lmidg98, jmid, fre, lbotznd, nas real ux12, uy12, uz12, ux13, uy13, uz13, unx, uny, unz, un real vx12, vy12, vz12, vx13, vy13, vz13, vnx, vny, vnz, vn real costh1, theta1, costh2, theta2, costh3, theta3 real x1, y1, z1, unc, v12, v13, d, co(3,mna) real lx, ly, lz, nrc, ql logical lprhb, flag, lprdu, doubleprint, flex character molcasin*(*), line*80 c ******************************************************************* nrc=real(normcharge) flag= .true. nfim= 1 do while (molcasin(nfim:nfim) .ne. ' ') nfim= nfim + 1 enddo nfim= nfim - 1 c escreve o output nos arquivos: *.xyz, *.or, *.dst c escreve primeiro o numero de atomos no arquivo *.xyz c write(6,*)'Solute:',nap,nmp if (nch .eq. 0) then write(uxyz, *) nap write(uall, *) nap else i= nmp+1 do while (i .le. mnm) k= hist(i) if (k .ne. 0) then nd = tim(k) nas = na(nd)-1 nap = nap + nas c write(6,*)'Achei cargas tipo:',nd,i,nas,nap,nmp+nch if (i .eq. (nmp+nch)) then i= mnm+1 else i= i + 1 endif else i= i + 1 endif enddo write(uxyz, *) nap write(uall, *) nap endif c escreve a linha de comentario no *.xyz write(uxyz,'('' Configuration number : '',i7,''L = '',3f9.4)') : iconf, lx, ly, lz write(uall,'('' Configuration number : '',i7,''L = '',3f9.4)') : iconf, lx, ly, lz if (iout .eq. 1) then C MOLCAS open(unit=umol, file= molcasin(:nfim), status='old') k= 1 j= tim(k) do while (flag) read(umol,'(a)',end=888) line j= index(line,'#orderch') if (j .ne. 0) then read(line(9:),*) nata ta= 0 do l= 2, na(j) if (ano(l,j) .eq. nata) then ta= ta + 1 if (flex) then co(1,l)= rx(k,l) co(2,l)= ry(k,l) co(3,l)= rz(k,l) else co(1,l)= rx1(l) co(2,l)= ry1(l) co(3,l)= rz1(l) endif if (ta .lt. 10) then write(uord,'(a,i1,3f15.8)') aicon(ano(l,j)),ta, : (co(il,l)/bohr,il=1,3) elseif (ta .lt. 100) then write(uord,'(a,i2,3f15.8)') aicon(ano(l,j)),ta, : (co(il,l)/bohr,il=1,3) else write(uord,'(a,i3,3f15.8)') aicon(ano(l,j)),ta, : (co(il,l)/bohr,il=1,3) endif endif enddo else write(uord,'(a)') line endif enddo 888 continue close(unit=umol) endif c escreve as coordenadas da primeira molecula k= 1 j= tim(k) do l= 2, na(j) if ((lprdu) .or. (ano(l,j) .ne. 104)) then if (flex) then co(1,l)= rx(k,l) co(2,l)= ry(k,l) co(3,l)= rz(k,l) else co(1,l)= rx1(l) co(2,l)= ry1(l) co(3,l)= rz1(l) endif if (iout .eq. 3) then C GAUSSIAN if (j.ne.bq) then write(uord,'(a,3f15.6)')aicon(ano(l,j)),(co(il,l),il=1,3) write(ug98,'(a,3f15.6)')aicon(ano(l,j)),(co(il,l),il=1,3) if ((confprint .eq. 1).and.(nmp.eq.1).and.(nch.gt.0) : .and.(fre.eq.1)) : write(ugpc,'(a,3f15.6)')aicon(ano(l,j)),(co(il,l),il=1,3) else c imprime com Bq write(uord,'(a,3f15.6)')aicon(105),(co(il,l),il=1,3) write(ug98,'(a,3f15.6)')aicon(105),(co(il,l),il=1,3) if ((confprint .eq. 1).and.(nmp.eq.1).and.(nch.gt.0) : .and.(fre.eq.1)) : write(ugpc,'(a,3f15.6)')aicon(105),(co(il,l),il=1,3) endif elseif (iout .ne. 1) then C ZINDO write(uord,'(3f10.5,i5)') (co(il,l),il=1,3), ano(l,j) if ((confprint .eq. 1).and.(nmp.eq.1).and.(nch.gt.0) : .and.(fre.eq.1)) : write(uzpc,'(3f10.5,i5)') (co(il,l),il=1,3), ano(l,j) endif write(uxyz,'(2x,a,3f15.5)')aicon(ano(l,j)),(co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)')aicon(ano(l,j)),(co(il,l),il=1,3) write(6,*) j, l, aicon(ano(l,j)), (co(il,l),il=1,3), lx,ly,lz if (ang .eq. 1) then l1= at1(j) + 1 l2= at2(j) + 1 l3= at3(j) + 1 x1 = rx(k,l1) y1 = ry(k,l1) z1 = rz(k,l1) ux12 = rx(k,l2) - x1 uy12 = ry(k,l2) - y1 uz12 = rz(k,l2) - z1 ux13 = rx(k,l3) - x1 uy13 = ry(k,l3) - y1 uz13 = rz(k,l3) - z1 unx = uy12*uz13 - uy13*uz12 uny = uz12*ux13 - uz13*ux12 unz = ux12*uy13 - ux13*uy12 unc = unx*x1 + uny*y1 + unz*z1 un = sqrt(unx*unx + uny*uny + unz*unz) endif endif enddo if (k .ne. 1) write(udst,'(2i5,3(2x,f9.5),2i4,2x,f9.6,2i4)') : iconf, k, sqrt(ri(k)),enrg(1,k), rmin1(k), timin1(k), : tjmin1(k), rmin2(k), timin2(k), tjmin2(k) ii= 1 c escreve as coordenadas das molecula ligadas por HB if (lprhb) then if (mathb .gt. 0) then do 100 i= 1, mathb k= athb(i) doubleprint= .false. if (i .gt. 1) then do i0= 1, i-1 if (k .eq. athb(i0)) doubleprint= .true. enddo endif j= tim(k) if ((ang .eq. 1) .and. (na(1) .gt. 3)) then l1= at1(j) + 1 l2= at2(j) + 1 l3= at3(j) + 1 x1 = rx(k,l1) y1 = ry(k,l1) z1 = rz(k,l1) d = unx*x1 + uny*y1 + unz*z1 if (d .lt. unc) then unx= -unx uny= -uny unz= -unz endif vx12 = rx(k,l2) - x1 vy12 = ry(k,l2) - y1 vz12 = rz(k,l2) - z1 vx13 = rx(k,l3) - x1 vy13 = ry(k,l3) - y1 vz13 = rz(k,l3) - z1 v12 = sqrt(vx12*vx12 + vy12*vy12 + vz12*vz12) v13 = sqrt(vx13*vx13 + vy13*vy13 + vz13*vz13) vnx = vy12*vz13 - vy13*vz12 vny = vz12*vx13 - vz13*vx12 vnz = vx12*vy13 - vx13*vy12 vn = sqrt(vnx*vnx + vny*vny + vnz*vnz) costh1= (unx*vnx + uny*vny + unz*vnz) / un / vn theta1= acos(costh1)* 180.0/pi costh2= (unx*vx12 + uny*vy12 + unz*vz12) / un / v12 theta2= acos(costh2)* 180.0/pi costh3= (unx*vx13 + uny*vy13 + unz*vz13) / un / v13 theta3= acos(costh3)* 180.0/pi if (k .ne. 1) write(udst,'(2i5,3(2x,f9.5),2i4,2x, : f9.6,2i4,3(3x,f7.2))') : iconf, k, sqrt(ri(k)),enrg(1,k), rmin1(k), timin1(k), : tjmin1(k),rmin2(k), timin2(k), tjmin2(k), theta1, : theta2, theta3 else if (k .ne. 1) write(udst,'(2i5,3(2x,f9.5),2i4,2x, : f9.6,2i4)') iconf, k, : sqrt(ri(k)),enrg(1,k), rmin1(k), timin1(k), tjmin1(k), : rmin2(k), timin2(k), tjmin2(k) endif if ((iout .eq. 1).and.(.not. doubleprint)) then C MOLCAS open(unit=umol, file= molcasin(:nfim), status='old') do while (flag) read(umol,'(a)',end=999) line j= index(line,'#orderch') if (j .ne. 0) then read(line(9:),*) nata do l= 2, na(j) if (ano(l,j) .eq. nata) then ta= ta + 1 co(1,l)= rx(k,l) co(2,l)= ry(k,l) co(3,l)= rz(k,l) if (ta .lt. 10) then write(uord,'(a,i1,3f15.8)') aicon(ano(l,j)), : ta,(co(il,l)/bohr,il=1,3) elseif (ta .lt. 100) then write(uord,'(a,i2,3f15.8)') aicon(ano(l,j)), : ta,(co(il,l)/bohr,il=1,3) else write(uord,'(a,i3,3f15.8)') aicon(ano(l,j)), : ta,(co(il,l)/bohr,il=1,3) endif endif enddo endif enddo 999 continue close(unit=umol) endif do 99 l= 2, na(j) if ((lprdu) .or. (ano(l,j) .ne. 104)) then co(1,l)= rx(k,l) co(2,l)= ry(k,l) co(3,l)= rz(k,l) if ((iout .eq. 3).and.(.not. doubleprint)) then C GAUSSIAN if (j.ne.bq) then write(uord,'(a,3f15.6)') aicon(ano(l,j)), : (co(il,l),il=1,3) write(ug98,'(a,3f15.6)') aicon(ano(l,j)), : (co(il,l),il=1,3) else c imprime com Bq write(uord,'(a,3f15.6)') aicon(105), : (co(il,l),il=1,3) write(ug98,'(a,3f15.6)') aicon(105), : (co(il,l),il=1,3) endif elseif ((iout .ne. 1).and.(.not. doubleprint)) then C ZINDO write(uord,'(3f10.5,i5)') (co(il,l),il=1,3),ano(l,j) endif write(uxyz,'(2x,a,3f15.5)') aicon(ano(l,j)), : (co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)') aicon(ano(l,j)), : (co(il,l),il=1,3) endif 99 continue ii= ii + 1 100 continue c escreve as coordenadas das cargas pontuais logo apos as c molecula ligadas por HB if (nch .gt. 0) then if (iout .eq. 1) then c MOLCAS write(uord,'(''Xfield'')') write(uord,*) nch*(na(2)-1) i = 2 jj= 1 do while (jj .le. nch) k = hist(i) do while (k .eq. 0) i= i+1 k= hist(i) enddo do while (hbflag(k)) i= i+1 k= hist(i) enddo j = tim(k) i = i + 1 jj= jj + 1 do l= 2, na(j) co(1,l)= rx(k,l) co(2,l)= ry(k,l) co(3,l)= rz(k,l) ql = z(l,j)/e if (abs(ql) .gt. 1E-5) then write(uord,'(3f15.8,x,f15.8,x,a)') : (co(il,l)/bohr,il=1,3),ql,' 0.000 0.000 0.000' endif write(uxyz,'(2x,a,3f15.5)')'Xx',(co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)')'Xx',(co(il,l),il=1,3) enddo enddo write(uord,'(''End of Input'')') elseif (iout .eq. 3) then c GAUSSIAN write(uord,*) write(ug98,*) c escreve o arquivo midfile logo apos as molecula ligadas c por HB e antes das cargas pontuais if (lmidg98 .gt. 0) then do jmid= 1, lmidg98 write(uord,'(a)') midg98(jmid) write(ug98,'(a)') midg98(jmid) enddo c write(ug98,*) endif i = 2 jj= 1 do while (jj .le. nch) k = hist(i) do while (k .eq. 0) i= i+1 k= hist(i) enddo do while (hbflag(k)) i= i+1 k= hist(i) enddo j = tim(k) i = i + 1 jj= jj + 1 do l= 2, na(j) co(1,l)= rx(k,l) co(2,l)= ry(k,l) co(3,l)= rz(k,l) ql = z(l,j)/e if (abs(ql) .gt. 1E-6) then write(uord,'(4f15.6)')(co(il,l),il=1,3),ql write(ug98,'(4f15.6)')(co(il,l),il=1,3),ql endif write(uxyz,'(2x,a,3f15.5)')'Xx',(co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)')'Xx',(co(il,l),il=1,3) enddo enddo else c ZINDO write(uord,'('' $END '')') write(uord,*) write(uord,'('' $PTCHGI '')') write(uord,'('' These are the point charges '')') write(uord,'('' 0.000000 0.000000 '')') i = 2 jj= 1 do while (jj .le. nch) k = hist(i) do while (k .eq. 0) i= i+1 k= hist(i) enddo do while (hbflag(k)) i= i+1 k= hist(i) enddo j = tim(k) i = i + 1 jj= jj + 1 do l= 2, na(j) co(1,l)= rx(k,l) co(2,l)= ry(k,l) co(3,l)= rz(k,l) ql = z(l,j)/e if (abs(ql) .gt. 1E-5) then write(uord,'(3f10.5,i5,f10.5)') : (co(il,l),il=1,3),0, ql endif write(uxyz,'(2x,a,3f15.5)')'Xx',(co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)')'Xx',(co(il,l),il=1,3) enddo enddo endif endif endif ii= ii - 1 c escreve as coordenadas demais molecula caso molprint > 1 ou =1 else if ((confprint.eq.1).and.(nmp.eq.1).and.(nch.gt.0) : .and.(fre.eq.1)) then if (iout .eq. 3) then c GAUSSIAN write(ugpc,*) c escreve o arquivo midfile no arquivo de configuracao c media antes das cargas pontuais if (lmidg98 .gt. 0) then do jmid= 1, lmidg98 write(ugpc,'(a)') midg98(jmid) enddo write(ugpc,*) endif elseif (iout .eq. 2) then c ZINDO write(uzpc,'('' $END '')') write(uzpc,*) write(uzpc,'('' $PTCHGI '')') write(uzpc,'(''These are the point charges '')') write(uzpc,'('' 0.000000 0.000000 '')') endif endif do 200 i= 2, mnm k= hist(i) if (k .ne. 0) then j= tim(k) if ((ang .eq. 1).and.(na(1) .gt. 3)) then l1= at1(j) + 1 l2= at2(j) + 1 l3= at3(j) + 1 x1 = rx(k,l1) y1 = ry(k,l1) z1 = rz(k,l1) d = unx*x1 + uny*y1 + unz*z1 if (d .lt. unc) then unx= -unx uny= -uny unz= -unz endif vx12 = rx(k,l2) - x1 vy12 = ry(k,l2) - y1 vz12 = rz(k,l2) - z1 vx13 = rx(k,l3) - x1 vy13 = ry(k,l3) - y1 vz13 = rz(k,l3) - z1 v12 = sqrt(vx12*vx12 + vy12*vy12 + vz12*vz12) v13 = sqrt(vx13*vx13 + vy13*vy13 + vz13*vz13) vnx = vy12*vz13 - vy13*vz12 vny = vz12*vx13 - vz13*vx12 vnz = vx12*vy13 - vx13*vy12 vn = sqrt(vnx*vnx + vny*vny + vnz*vnz) costh1= (unx*vnx + uny*vny + unz*vnz) / un / vn theta1= acos(costh1)* 180.0/pi costh2= (unx*vx12 + uny*vy12 + unz*vz12) / un / v12 theta2= acos(costh2)* 180.0/pi costh3= (unx*vx13 + uny*vy13 + unz*vz13) / un / v13 theta3= acos(costh3)* 180.0/pi if (k .ne. 1) write(udst,'(2i5,3(2x,f9.5),2i4,2x, : f9.6,2i4,3(3x,f7.2))') : iconf, k, sqrt(ri(k)),enrg(1,k), rmin1(k), timin1(k), : tjmin1(k), rmin2(k), timin2(k), tjmin2(k), theta1, : theta2, theta3 else if (k .ne. 1) write(udst,'(2i5,3(2x,f9.5),2i4,2x, : f9.6,2i4)') : iconf, k, sqrt(ri(k)),enrg(1,k), rmin1(k), timin1(k), : tjmin1(k), rmin2(k), timin2(k), tjmin2(k) endif ii= ii + 1 do 199 l= 2, na(j) co(1,l)= rx(k,l) co(2,l)= ry(k,l) co(3,l)= rz(k,l) if (iout .eq. 1) then C MOLCAS if (ii .le. nmp) then if ((lprdu) .or. (ano(l,j) .ne. 104)) : write(uord,'(a,x,3f15.8)')aicon(ano(l,j)), : (co(il,l)/bohr,il=1,3) else if (nch .gt. 0) then if ((ii .eq. (nmp+1)) .and. (l .eq. 2)) then write(uord,'(''Xfield'')') write(uord,*) nch*(na(2)-1) endif ql = z(l,j)/e if (abs(ql) .gt. 1E-8) then write(uord,'(3f15.8,x,f15.8,x,a)') : (co(il,l)/bohr,il=1,3), ql, ' 0.000 0.000 :0.000' endif write(uxyz,'(2x,a,3f15.5)') 'Xx', : (co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)') 'Xx', : (co(il,l),il=1,3) endif endif elseif (iout .eq. 3) then C GAUSSIAN if (ii .le. nmp) then if ((lprdu) .or. (ano(l,j) .ne. 104)) then if (j.ne.bq) then write(uord,'(a,3f15.6)') aicon(ano(l,j)), : (co(il,l),il=1,3) write(ug98,'(a,3f15.6)') aicon(ano(l,j)), : (co(il,l),il=1,3) else c imprime Bq write(uord,'(a,3f15.6)') aicon(105), : (co(il,l),il=1,3) write(ug98,'(a,3f15.6)') aicon(105), : (co(il,l),il=1,3) endif endif else if (nch .gt. 0) then if ((ii .eq. (nmp+1)) .and. (l .eq. 2)) then write(uord,*) write(ug98,*) c escreve o arquivo midfile antes das cargas pontuais if (lmidg98 .gt. 0) then do jmid= 1, lmidg98 write(uord,'(a)') midg98(jmid) write(ug98,'(a)') midg98(jmid) enddo c write(ug98,*) endif endif ql = z(l,j)/e if (abs(ql) .gt. 1E-6) then write(uord,'(4f15.6)')(co(il,l),il=1,3),ql write(ug98,'(4f15.6)')(co(il,l),il=1,3),ql if ((nmp.eq.1).and.(fre.eq.1)) : write(ugpc,'(4f15.6)')(co(il,l),il=1,3),ql/nrc endif write(uxyz,'(2x,a,3f15.5)') 'Xx', : (co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)') 'Xx', : (co(il,l),il=1,3) endif endif else C ZINDO if (ii .le. nmp) then if ((lprdu) .or. (ano(l,j) .ne. 104)) : write(uord,'(3f10.5,i5)') (co(il,l),il=1,3),ano(l,j) else if (nch .gt. 0) then if((ii .eq. (nmp+1)) .and. (l .eq. 2)) then write(uord,'('' $END '')') write(uord,*) write(uord,'('' $PTCHGI '')') write(uord,'(''These are the point charges '')') write(uord,'('' 0.000000 0.000000 '')') endif ql = z(l,j)/e if (abs(ql) .gt. 1E-5) then write(uord,'(3f10.5,i5,f10.5)') : (co(il,l),il=1,3),0, ql if ((nmp.eq.1).and.(fre.eq.1)) : write(uzpc,'(3f10.5,i5,f10.5)') : (co(il,l),il=1,3),0, ql/nrc endif write(uxyz,'(2x,a,3f15.5)') 'Xx', : (co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)') 'Xx', : (co(il,l),il=1,3) endif endif endif if (ii .le. nmp) then if ((lprdu) .or. (ano(l,j) .ne. 104)) then write(uxyz,'(2x,a,3f15.5)') aicon(ano(l,j)), : (co(il,l),il=1,3) write(uall,'(2x,a,3f15.5)') aicon(ano(l,j)), : (co(il,l),il=1,3) endif endif 199 continue endif if (ii .eq. (nmp+nch)) goto 900 200 continue endif 900 continue if (iout .eq. 1) then C MOLCAS write(uord,'(''End of Input'')') elseif (iout .eq. 3) then C GAUSSIAN c escreve o arquivo midfile logo apos as N moleculas c se numa carga pontual for escrita if ((lmidg98 .gt. 0).and.(nch.eq.0)) then write(uord,*) write(ug98,*) do jmid= 1, lmidg98 write(uord,'(a)') midg98(jmid) write(ug98,'(a)') midg98(jmid) enddo c write(ug98,*) endif write(ug98,'('' '')') write(ug98,'(''--link1--'')') else c ZINDO if (lbotznd .eq. 0) then write(uord,'('' $END '')') write(uord,*) endif endif return end