!********************************************************** ! SUBROUTINE : khydration ! ! Purpose : return the list of hydrate (1 water molecule added) ! that can be made from the formula given as input) ! ! INPUT : ! - chem : formula of the input species ! - yield : ratio of the formula. Species may have multi hydration. ! The yield provide the ratio with respect to the fully ! non hydrated molecule. ! ! Output : ! - ndat : number of distinct molecule (hydrate) that can be made ! from chem (again 1 water molecule added only !) ! - chemdat(i) : formula of the i'th hydrate ! - ydat(i) : hydration constant of the i'th hydrate (with respect ! to the fully non hydrated molecule !) ! ! ! INTERNAL ! - nabcde(k) : number of distinct pathways that end up at a ! position k relative to top (e.g. nabcde(3) gives ! the number of distinct pathways finishing in a ! beta position relative to top ! - tabcde(k,i,j) : give the pathways (nodes j), for the track number i ! to reach the position k (k=2 is beta position ...). ! For example, tabcde(4,1,j) give the first track to ! reach a gamma position (node given by ! tabcde(4,1,4), using the track given by ! tabcde(4,1,*) ! - mapfun(a,b,c) : provide the number of function of type 'c' ! at position (node) 'a'. index 'b' if for node ! type with 1=aliphatic, 2=cd and 3=aromatic ! for example, the molecule CH2(OH)CdH=CdHCHO ! should have non zero values at the positions : ! mapfun(1,1,1)=1 ! mapfun(4,2,9)=1 ! - funflg(a) : get the number of functional group at ! node a. For the example above, non-zero ! values are found at position 1 and 4, where ! it is set to 1. ! - nodetype : table of character for type node ! 'y' is for carbonyl ! 'r' is for aromatic ! 'o' is for -O- node ! 'd' is for Cd ! 'n' is for the others (i.e. normal) ! ! * ! table of functions - Index of functionalities (in mapfun) ! ----------------------------------------------------------------- ! 1= -OH ; 2= -NO2 ; 3= -ONO2 ; 4= -OOH ; 5= -F ! 6= -Cl ; 7= -Br ; 8= -I ; 9= -CHO ; 10= -CO- ! 11= -COOH ; 12= -CO(OOH) ; 13= -PAN ; 14= -O- ; 15= R-COO-R ! 16 = HCO-O-R; 17= -CO(F) ; 18= -CO(Cl) ; 19= -CO(Br) ; 20= -CO(I) ! This list is used in alisig which provide sigma for a given group ! ---------------------------------------------------------------------- !********************************************************** SUBROUTINE khydration(chem,yield,ndat,chemdat,ydat) IMPLICIT NONE INCLUDE 'general.h' ! input CHARACTER(LEN=lfo) chem REAL yield ! output INTEGER ndat CHARACTER(LEN=lfo) chemdat(mhyd) REAL ydat(mhyd) ! local CHARACTER(LEN=lfo) tchemdat(mhyd) REAL tydat(mhyd) INTEGER tndat INTEGER nc,i,j,k,l INTEGER bond(mca,mca), dbflg, nring, node CHARACTER(LEN=lgr) group(mca) CHARACTER*1 nodetype(mca) REAL alifun(20),cdfun(20),arofun(20) REAL mapfun(mca,3,20) INTEGER funflg(mca) INTEGER tabester(4,2) ! 1= -O- side, 2= CO side INTEGER ngrp INTEGER ierr CHARACTER(LEN=lfo) tempkc CHARACTER(LEN=lgr) tgroup(mca), pold, pnew INTEGER tbond(mca,mca) INTEGER rjg(mri,2) ! ring-join group pairs INTEGER nabcde(9), tabcde(9,mco,mca) REAL alisig(20),arometa(20),aropara(20),arorto(20) REAL sigester REAL sigma, hsigma REAL khyd, sumy INTEGER numy INTEGER nfcd,nfcr ! ----------- ! initialize ! ----------- nfcr=0 nfcd=0 ndat=0 tndat=0 DO i=1,mhyd chemdat(i)=' ' ydat(i)=0. tchemdat(i)=' ' tydat(i)=0. ENDDO ! build the group and bond matrix for chem nc=index(chem,' ') CALL grbond(chem,nc,group,bond,dbflg,nring) DO i=1,mca bond(i,i)=0 ENDDO CALL rjgrm(nring,group,rjg) ! rm ring index and get ring closing nodes ! make copy DO i=1,mca tgroup(i)=group(i) DO j=1,mca tbond(i,j)=bond(i,j) ENDDO ENDDO ! count the number of group in chem and get the map of the functional ! group (see index number related the organic function) node=0 DO i=1,mca IF (group(i)(1:1).NE.' ') node = node + 1 ENDDO CALL chemmap(chem,node,group,bond,ngrp,nodetype, & alifun,cdfun,arofun,mapfun,funflg, & tabester,nfcd,nfcr,ierr) IF (ierr.ne.0) STOP 'stop in khydration' IF (ierr.ne.0) RETURN ! --------------------------------- ! hydration of aliphatic aldehyde ! --------------------------------- IF ((alifun(9).NE.0).OR. & (cdfun(9).NE.0) )THEN ! find aldehyde position DO 10 i=1,node IF ( (mapfun(i,1,9).eq.1.).OR. & (mapfun(i,2,9).eq.1.) )THEN sigma=0. CALL abcde_map(bond,i,node,nabcde,tabcde) ! get the neighboors. ! check that the aldehyde is not an hidden anhydride (-CO-O-CO-) DO k=1,nabcde(2) IF (nodetype(tabcde(2,k,2)).eq.'o') GOTO 10 ENDDO ! Remove the path that would lead to count the ester twice. ! Performed by "killing" the "second" node of the functionality. CALL kill_ester(tabester,nabcde,tabcde) ! get taft sigma (look for groups up to delta position) CALL get_tsig(bond,group,5,nabcde,tabcde,mapfun, & funflg,nodetype,sigma) ! compute the hydration constant for aliphatic aldehyde khyd=0.0818 + 1.2663*sigma ! If the species is a 'Cd' aldehyde, then add the correction applied ! for aromatic (no data available) IF (mapfun(i,2,9).eq.1.) khyd=khyd-1.5817 ! make the hydrate, rebuild, rename and reset tgroup & tbond pold = 'CHO ' pnew = 'CH(OH)(OH) ' CALL swap(group(i),pold,tgroup(i),pnew) CALL rebond(tbond,tgroup,tempkc,nring) ! print*,tempkc CALL stdchm(tempkc) DO k=1,mca tgroup(k)=group(k) DO l=1,mca tbond(k,l)=bond(k,l) ENDDO ENDDO ! store the data tndat=tndat+1 IF (tndat.gt.mhyd) STOP 'in khydration, exceed max hydrate' tydat(tndat)=khyd+yield tchemdat(tndat)=tempkc ENDIF 10 CONTINUE ENDIF ! end of aliphatic aldehyde ! --------------------------------- ! hydration of aliphatic ketone ! --------------------------------- IF (alifun(10).NE.0) THEN ! find ketone position DO 20 i=1,node IF (mapfun(i,1,10).eq.1.) THEN sigma=0. ! get the neighboors. CALL abcde_map(bond,i,node,nabcde,tabcde) ! check that the ketone is not an hidden anhydride (-CO-O-CO-) DO k=1,nabcde(2) IF (nodetype(tabcde(2,k,2)).eq.'o') GOTO 20 ENDDO ! Remove the path that would lead to count the ester twice. ! Performed by "killing" the "second" node of the functionality. CALL kill_ester(tabester,nabcde,tabcde) ! get taft sigma (look for groups up to delta position) CALL get_tsig(bond,group,5,nabcde,tabcde,mapfun, & funflg,nodetype,sigma) ! compute the hydration constant for aliphatic ketone khyd=0.0818 + 1.2663*sigma - 2.5039 ! make the hydrate, rebuild, rename and reset tgroup & tbond pold = 'CO ' pnew = 'C(OH)(OH) ' CALL swap(group(i),pold,tgroup(i),pnew) CALL rebond(tbond,tgroup,tempkc,nring) CALL stdchm(tempkc) DO k=1,mca tgroup(k)=group(k) DO l=1,mca tbond(k,l)=bond(k,l) ENDDO ENDDO ! store the data tndat=tndat+1 IF (tndat.gt.mhyd) STOP 'in khydration, exceed max hydrate' tydat(tndat)=khyd+yield tchemdat(tndat)=tempkc ENDIF 20 CONTINUE ENDIF ! end of aliphatic ketone ! --------------------------------- ! hydration of aromatic aldehyde ! --------------------------------- IF (arofun(9).NE.0) THEN ! find aldehyde position DO i=1,node IF (mapfun(i,3,9).eq.1.) THEN sigma=0. CALL abcde_map(bond,i,node,nabcde,tabcde) ! get the neighboors. ! Remove the path that would lead to count the ester twice. ! Performed by "killing" the "second" node of the functionality. CALL kill_ester(tabester,nabcde,tabcde) ! get hammet sigma (look for groups up to delta position) CALL get_hsig(bond,group,nabcde,tabcde,mapfun, & funflg,nodetype,hsigma) ! compute the hydration constant for alihpatic aldehyde khyd=0.0818 + 0.4969*hsigma - 1.5817 ! make the hydrate, rebuild, rename and reset tgroup & tbond pold = 'CHO ' pnew = 'CH(OH)(OH) ' CALL swap(group(i),pold,tgroup(i),pnew) CALL rebond(tbond,tgroup,tempkc,nring) CALL stdchm(tempkc) DO k=1,mca tgroup(k)=group(k) DO l=1,mca tbond(k,l)=bond(k,l) ENDDO ENDDO ! store the data tndat=tndat+1 IF (tndat.gt.mhyd) STOP 'in khydration, exceed max hydrate' tydat(tndat)=khyd+yield tchemdat(tndat)=tempkc ENDIF ENDDO ENDIF ! end of aromatic aldehyde ! ---------------------------------------------- ! collapse identical product ! ---------------------------------------------- ! NOTE : must add khyd, not the log of it DO i=1,tndat numy=1 sumy=10**(tydat(i)) DO j=i+1,tndat IF (tchemdat(i).eq.tchemdat(j)) THEN numy=numy+1 sumy=sumy+10**(tydat(j)) tchemdat(j)=' ' ENDIF ENDDO IF (numy.gt.1 .and. sumy .gt. 0) THEN tydat(i)=log10(sumy) ENDIF ENDDO DO i=1,tndat IF (tchemdat(i)(1:1).ne.' ') THEN ndat=ndat+1 chemdat(ndat)=tchemdat(i) ydat(ndat)=tydat(i) ENDIF ENDDO END