**************************************************************************** * * * Estimate the boiling points with the SIMPOL method (Pankow/Asher 2008) * * https://doi.org/10.5194/acp-8-2773-2008 * * **************************************************************************** SUBROUTINE simpol1(chem,group,bond,nring,simpgroup) IMPLICIT NONE INCLUDE 'general.h' * input: CHARACTER(LEN=lfo) chem CHARACTER(LEN=lgr) group(mca) INTEGER bond(mca,mca),nring REAL temp * output REAL simpgroup(0:30) * internal CHARACTER(LEN=lgr) neigh(mca,4) REAL bk(0:30,4) INTEGER i,j,k INTEGER begrg,endrg INTEGER ring(mca) ! rg index for nodes, current ring (0=no, 1=yes) INTEGER indexrg(mca) ! rg index for nodes, any ring (0=no, 1=yes) INTEGER rngflg ! 0 = 'no ring', 1 = 'yes ring' INTEGER rjg(mri,2) ! ring-join group pairs INTEGER nca ! number of nodes INTEGER nc, ibeg, eflg INTEGER nbnei(mca),nei_ind(mca,4),nbatom REAL summ INTEGER ncd,conjug INTEGER ncdcase, cdtable(4),tcdtable(4), cdsub(4) INTEGER cdcarbo(4,2),cdeth(4,2),lring(2) INTEGER rxnflg REAL logPvap,dH *---------------------------------------- * initialization *---------------------------------------- * count the number of nodes nca=0 ring=0 nbnei=0 DO i=1,mca IF (group(i)(1:1).NE.' ') nca=nca+1 ENDDO neigh=' ' simpgroup=0 bk=0 rxnflg = 0 * IF RINGS EXIST, THEN FIND THE NODES THAT BELONG TO THE RINGS * -------------------------------------------------------------- IF (nring.gt.0) THEN * remove ring-join characters from groups and find the nodes that close rings CALL rjgrm(nring,group,rjg) * find the nodes that belong to a ring DO i=1,nring begrg=rjg(i,1) endrg=rjg(i,2) CALL findring(begrg,endrg,nca,bond,rngflg,ring) ENDDO ENDIF * find neighbours DO i=1,nca DO j=1,nca IF (bond(i,j).NE.0) THEN nbnei(i)=nbnei(i)+1 !number of neighbours neigh(i,nbnei(i))=group(j) !groups of neighbours nei_ind(i,nbnei(i))=j ENDIF ENDDO ENDDO *---------------------------------------- * sum of contributions for each group *---------------------------------------- simpgroup(0)=1 simpgroup(1)=nca DO i=1,nca IF (group(i).EQ.'-O-') simpgroup(1)=simpgroup(1)-1 ENDDO * test aromatic ring IF (INDEX(chem,'c1').NE.0) THEN simpgroup(3)=nring ELSE simpgroup(4)=nring ENDIF DO i=1,nca IF (group(i)(1:2).EQ.'Cd') THEN simpgroup(5) = simpgroup(5) + 0.5 CALL cdcase_gt2(chem,bond,group,rxnflg,nring, & ncd,ncdcase,conjug, & cdtable,tcdtable,cdsub,cdcarbo,cdeth) !IF (ncdcase.EQ.2) THEN ! Simpgroup(6) C=C-C=O in aliphatic ring IF (ncdcase.EQ.2.AND. & ring(i).NE.0.AND.INDEX(group(i),"c").EQ.0) THEN simpgroup(6) = simpgroup(6) + 0.5 ENDIF ENDIF ! If only ONE of the Cd carbons is on a ring, group(6) does not apply. simpgroup(6) = INT(simpgroup(6)) IF ((INDEX(group(i),'(OH)').NE.0).AND. & (group(i).NE.'CO(OH)')) THEN IF (group(i)(1:1).NE.'c') THEN simpgroup(7)=simpgroup(7)+1 ELSE simpgroup(17)=simpgroup(17)+1 ENDIF ENDIF IF ((INDEX(group(i),'(OOH)').NE.0).AND. & (group(i).NE.'CO(OOH)')) THEN simpgroup(27)=simpgroup(27)+1 ENDIF IF (group(i).EQ.'CHO') THEN simpgroup(8)=simpgroup(8)+1 ENDIF IF (group(i).EQ.'CO') THEN IF ((neigh(i,1).EQ.'-O-').AND.(neigh(i,2).EQ.'-O-')) THEN simpgroup(12)=simpgroup(12)+1 GOTO 10 ELSE IF ((neigh(i,1).EQ.'-O-').OR.(neigh(i,2).EQ.'-O-')) THEN simpgroup(11)=simpgroup(11)+0.5 GOTO 10 ENDIF simpgroup(9)=simpgroup(9)+1 ENDIF 10 CONTINUE IF (group(i).EQ.'-O-') THEN IF ((neigh(i,1).EQ.'CO').AND.(neigh(i,2).EQ.'CO')) THEN simpgroup(9)=simpgroup(9)+1 GOTO 20 ELSE IF ((neigh(i,1).EQ.'CO').OR.(neigh(i,2).EQ.'CO')) THEN simpgroup(11)=simpgroup(11)+0.5 GOTO 20 ENDIF IF ((neigh(i,1)(1:1).NE.'c').AND.(neigh(i,2)(1:1).NE.'c') & .AND.(ring(i).EQ.0)) THEN simpgroup(12)=simpgroup(12)+1 ELSE IF (ring(i).NE.0) THEN simpgroup(13)=simpgroup(13)+1 ELSE simpgroup(14)=simpgroup(14)+1 ENDIF ENDIF 20 CONTINUE IF (group(i).EQ.'CO(OH)') THEN simpgroup(10)=simpgroup(10)+1 ENDIF IF (group(i).EQ.'CO(OOH)') THEN simpgroup(28)=simpgroup(28)+1 ENDIF IF (group(i).EQ.'CO(OONO2)') THEN simpgroup(25)=simpgroup(25)+1 ENDIF IF (INDEX(group(i),'(ONO2)').NE.0) THEN simpgroup(15)=simpgroup(15)+1 ENDIF IF (INDEX(group(i),'(NO2)').NE.0) THEN simpgroup(16)=simpgroup(16)+1 IF (((neigh(i,1).EQ.'CO').OR.(neigh(i,2).EQ.'CO')) & .AND.(simpgroup(11).NE.0)) THEN simpgroup(30)=simpgroup(30)+1 simpgroup(11)=simpgroup(11)-1 ENDIF ENDIF ENDDO DO i=1,nca IF ((simpgroup(16).NE.0).AND.(simpgroup(17).NE.0)) THEN simpgroup(29)=simpgroup(17) simpgroup(17)=0. ENDIF ENDDO ********************************************************** * check errors * ********************************************************** DO i=0,30 IF (simpgroup(i).LT.0) THEN WRITE(6,*) '--error in simpol1.f --' WRITE(6,*) 'simpgroup negatif - see fort.43' WRITE(6,*) 'simpgroup(',i,')=',simpgroup(i) WRITE(43,*) 'chem=',chem(1:50) DO j=0,30 IF (simpgroup(j).NE.0) THEN WRITE(43,*) 'simpgroup(',j,')=',simpgroup(j) ENDIF ENDDO STOP ENDIF summ=simpgroup(i)*2 IF (mod(summ,2.).NE.0) THEN WRITE(6,*) '--error in simpol1.f --' WRITE(6,*) 'number on simpgroup not an integer' WRITE(6,*) 'simpgroup(',i,')=',simpgroup(i) WRITE(43,*) 'chem=',chem(1:50) DO j=0,30 IF (simpgroup(j).NE.0) THEN WRITE(43,*) 'simpgroup(',j,')=',simpgroup(j) ENDIF ENDDO STOP ENDIF ENDDO *---------------------------------------- * reading data *---------------------------------------- c! OPEN(10,file='../DATA/simpol.dat',form='FORMATTED', status='OLD') c OPEN(10,file='../DATA/simpol.dat',status='OLD') c DO i=0,30 c READ(10,'(f10.4,2x,f9.6,2x,f11.8,2x,f11.8))') c & bk(i,1),bk(i,2),bk(i,3),bk(i,4) c ENDDO c CLOSE(10) *---------------------------------------- * writing output *---------------------------------------- logPvap=0 dH=0 temp=298. DO i=0,30 IF (simpgroup(i).NE.0) THEN logPvap = logPvap + simpgroup(i)* & ((bk(i,1)/temp)+bk(i,2)+(bk(i,3)*temp)+(bk(i,4)*alog(temp))) dH=dH + simpgroup(i)* & ( bk(i,1) - (bk(i,3)*temp*temp) - (bk(i,4)*temp)) ENDIF ENDDO dH= -2.303 * 8.314 * dH c DO i=0,30 c IF (simpgroup(i).NE.0) THEN cc WRITE(6,*) 'simpgroup(',i,')=',simpgroup(i) c ENDIF c ENDDO c WRITE(57,'(A50,2x,31(f3.0))') chem(1:50),(simpgroup(i),i=0,30) RETURN END