! #include "copyright.h" module genhpmf !+++++++++++++++++++++++++++++++++++++++++++ ! ! ----Module for HPMF called genhpmf. ! ! This module contains main function ehpmf, ! subroutines sort, init_hpmf,attach, khpmf !++++++++++++++++++++++++++++++++++++++++++++ #include "assert.h" #include "dprec.h" #include "memory.h" integer,private :: n integer,private,parameter :: maxatm = 25000 integer,private :: atomic(maxatm) _REAL_,private :: x(maxatm) _REAL_,private :: y(maxatm) _REAL_,private :: z(maxatm) integer,private,parameter :: maxval = 30 integer,private,parameter :: maxn13 = 3*maxval integer,private,parameter :: maxn14 = 9*maxval integer,private,parameter :: maxn15 = 27*maxval integer,private :: n12(maxatm) integer,private,dimension(maxval,maxatm) :: i12 integer,private :: n13(maxatm) integer,private,dimension(maxn13,maxatm) :: i13 integer,private :: n14(maxatm) integer,private,dimension(maxn14,maxatm) :: i14 integer,private :: n15(maxatm) integer,private,dimension(maxn15,maxatm) :: i15 _REAL_,private,parameter :: rcarbon=1.7d0 _REAL_,private,parameter :: rwater=1.4d0 _REAL_,private,parameter :: acsurf=120.7628d0 _REAL_,private,parameter :: safact=0.3516d0 _REAL_,private,parameter :: tslope=100.0d0 _REAL_,private,parameter :: toffset=6.0d0 _REAL_,private,parameter :: hpmfcut=11.0d0 _REAL_,private :: h1,h2,h3 _REAL_,private :: c1,c2,c3 _REAL_,private :: w1,w2,w3 _REAL_,private :: radian _REAL_,private :: pi _REAL_,private :: sqrtpi _REAL_,private :: logten _REAL_,private :: sqrttwo _REAL_,private :: twosix parameter(h1=-0.7308004860404441194d0) parameter(h2=0.2001645051578760659d0) parameter(h3=-0.0905499953418473502d0) parameter(c1=3.8167879266271396155d0) parameter(c2=5.4669162286016419472d0) parameter(c3=7.1167694861385353278d0) parameter(w1=1.6858993102248638341d0) parameter(w2=1.3906405621629980285d0) parameter(w3=1.5741657341338335385d0) parameter (radian=57.29577951308232088d0) parameter (pi=3.141592653589793238d0) parameter (sqrtpi=1.772453850905516027d0) parameter (logten=2.302585092994045684d0) parameter (sqrttwo=1.414213562373095049d0) parameter (twosix=1.122462048309372981d0) integer,private :: npmf integer,private :: ipmf(maxatm) _REAL_,private :: rpmf(maxatm) _REAL_,private :: acsa(maxatm) contains !+++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! --- assign the variables of AMBER to HPMF in TINKER ! !++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine init_hpmf(xx,ix,ih,natom,nbona) implicit none integer natom,nbona _REAL_ xx(*) integer ix(*) character(len=4) ih(*) character(len=2) atype integer i,j integer jj,kk n = natom do i=1,n x(i) = xx(3*i-2) y(i) = xx(3*i-1) z(i) = xx(3*i) read(ih(m06+i-1),*) atype if(atype(1:1)=='H' .or. atype(1:1)=='h') atomic(i) = 1 if(atype(1:1)=='B' .or. atype(1:1)=='b') atomic(i) = 5 if(atype(1:1)=='C' .or. atype(1:1)=='c') atomic(i) = 6 if(atype(1:1)=='N' .or. atype(1:1)=='n') atomic(i) = 7 if(atype(1:1)=='O' .or. atype(1:1)=='o') atomic(i) = 8 if(atype(1:1)=='F' .or. atype(1:1)=='f') atomic(i) = 9 if(atype(1:1)=='P' .or. atype(1:1)=='p') atomic(i) = 15 if(atype(1:1)=='S' .or. atype(1:1)=='s') atomic(i) = 16 if(atype(1:1)=='I' .or. atype(1:1)=='i') atomic(i) = 53 if(atype == 'CT') atomic(i) = 6 if(atype == '1H' .or. atype == '2H' .or. atype == '3H') atomic(i) = 1 if(atype == 'HE' .or. atype == 'he') atomic(i) = 2 if(atype == 'NE' .or. atype == 'ne') atomic(i) = 10 if(atype == 'SI' .or. atype == 'si') atomic(i) = 14 if(atype == 'CL' .or. atype == 'cl') atomic(i) = 17 if(atype == 'AR' .or. atype == 'ar') atomic(i) = 18 if(atype == 'SE' .or. atype == 'se') atomic(i) = 34 if(atype == 'BR' .or. atype == 'br') atomic(i) = 35 if(atype == 'KR' .or. atype == 'kr') atomic(i) = 36 if(atype == 'XE' .or. atype == 'xe') atomic(i) = 54 end do do i=1,natom n12(i) = 0 end do do i= 1, nbona jj = ix(iiba+3)/3 kk = ix(ijba+3)/3 n12(jj) = n12(jj) + 1 i12(n12(jj),jj) = kk n12(kk) = n12(kk) + 1 i12(n12(kk),kk) = jj end do do i=1,natom call sort (n12(i),i12(1,i)) end do return end subroutine init_hpmf !+++++++++++++++++++++++++++++ ! ! "khpmf" initializes parameters needed for the hydrophobic ! potential of mean force nonpolar implicit solvation model1 ! !+++++++++++++++++++++++++++++ subroutine khpmf implicit none integer i,j,k integer nh,atmnum logical keep ! get carbons for PMF and set surface area screening values npmf = 0 do i = 1, n if (atomic(i) .eq. 6) then keep = .true. nh = 0 if (n12(i) .le. 2) keep = .false. do j = 1, n12(i) k = i12(j,i) if (atomic(k) .eq. 1) nh = nh + 1 if (n12(i).eq.3 .and. atomic(k).eq.8) keep = .false. end do if (keep) then npmf = npmf + 1 ipmf(npmf) = i acsa(i) = 1.0d0 if (n12(i).eq.3 .and. nh.eq.0) acsa(i) = 1.554d0 if (n12(i).eq.3 .and. nh.eq.1) acsa(i) = 1.073d0 if (n12(i).eq.4 .and. nh.eq.1) acsa(i) = 1.276d0 if (n12(i).eq.4 .and. nh.eq.2) acsa(i) = 1.045d0 if (n12(i).eq.4 .and. nh.eq.3) acsa(i) = 0.880d0 acsa(i) = acsa(i) * safact/acsurf end if end if end do ! assign HPMF atomic radii from traditional Bondi values do i = 1, n rpmf(i) = 2.0d0 atmnum = atomic(i) if (atmnum .eq. 0) rpmf(i) = 0.00d0 if (atmnum .eq. 1) rpmf(i) = 1.20d0 if (atmnum .eq. 2) rpmf(i) = 1.40d0 if (atmnum .eq. 5) rpmf(i) = 1.80d0 if (atmnum .eq. 6) rpmf(i) = 1.70d0 if (atmnum .eq. 7) rpmf(i) = 1.55d0 if (atmnum .eq. 8) rpmf(i) = 1.50d0 if (atmnum .eq. 9) rpmf(i) = 1.47d0 if (atmnum .eq. 10) rpmf(i) = 1.54d0 if (atmnum .eq. 14) rpmf(i) = 2.10d0 if (atmnum .eq. 15) rpmf(i) = 1.80d0 if (atmnum .eq. 16) rpmf(i) = 1.80d0 if (atmnum .eq. 17) rpmf(i) = 1.75d0 if (atmnum .eq. 18) rpmf(i) = 1.88d0 if (atmnum .eq. 34) rpmf(i) = 1.90d0 if (atmnum .eq. 35) rpmf(i) = 1.85d0 if (atmnum .eq. 36) rpmf(i) = 2.02d0 if (atmnum .eq. 53) rpmf(i) = 1.98d0 if (atmnum .eq. 54) rpmf(i) = 2.16d0 end do return end subroutine khpmf !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! "attach" generates lists of 1-3, 1-4 and 1-5 connectivities ! starting from the previously determined list of attached ! atoms (ie, 1-2 connectivity) ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine attach implicit none integer i,j,k,m integer jj,kk ! loop over all atoms finding all the 1-3 relationships; ! note "n12" and "i12" have already been setup elsewhere do i = 1, n n13(i) = 0 do j = 1, n12(i) jj = i12(j,i) do k = 1, n12(jj) kk = i12(k,jj) if (kk .eq. i) goto 10 do m = 1, n12(i) if (kk .eq. i12(m,i)) goto 10 end do n13(i) = n13(i) + 1 i13(n13(i),i) = kk 10 continue end do end do call sort (n13(i),i13(1,i)) end do ! loop over all atoms finding all the 1-4 relationships do i = 1, n n14(i) = 0 do j = 1, n13(i) jj = i13(j,i) do k = 1, n12(jj) kk = i12(k,jj) if (kk .eq. i) goto 30 do m = 1, n12(i) if (kk .eq. i12(m,i)) goto 30 end do do m = 1, n13(i) if (kk .eq. i13(m,i)) goto 30 end do n14(i) = n14(i) + 1 i14(n14(i),i) = kk 30 continue end do end do call sort (n14(i),i14(1,i)) end do ! loop over all atoms finding all the 1-5 relationships do i = 1, n n15(i) = 0 do j = 1, n14(i) jj = i14(j,i) do k = 1, n12(jj) kk = i12(k,jj) if (kk .eq. i) goto 50 do m = 1, n12(i) if (kk .eq. i12(m,i)) goto 50 end do do m = 1, n13(i) if (kk .eq. i13(m,i)) goto 50 end do do m = 1, n14(i) if (kk .eq. i14(m,i)) goto 50 end do n15(i) = n15(i) + 1 i15(n15(i),i) = kk 50 continue end do end do call sort (n15(i),i15(1,i)) end do return end subroutine attach !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! !"ehpmf" calculates the hydrophobic potential of mean force ! energy using a pairwise double loop ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ _REAL_ function ehpmf() implicit none integer i,j,k,m integer ii,jj,kk integer sschk integer omit(maxatm) _REAL_ xr,yr,zr,r,r2 _REAL_ e,ehp _REAL_ rsurf,pisurf _REAL_ hpmfcut2 _REAL_ saterm,sasa _REAL_ rbig,rsmall _REAL_ part,cutv _REAL_ e1,e2,e3,sum _REAL_ arg1,arg2,arg3 _REAL_ arg12,arg22,arg32 _REAL_ cutmtx(maxatm) call attach call khpmf ! zero out the hydrophobic potential of mean force energy ehp = 0.0d0 ! set some values needed during the HPMF calculation rsurf = rcarbon + 2.0d0*rwater pisurf = pi * (rcarbon+rwater) hpmfcut2 = hpmfcut * hpmfcut ! get the solvent accessible surface area for each atom do ii = 1, npmf i = ipmf(ii) saterm = acsa(i) sasa = 1.0d0 do k = 1, n if (i .ne. k) then xr = x(i) - x(k) yr = y(i) - y(k) zr = z(i) - z(k) r2 = xr*xr + yr*yr + zr*zr rbig = rpmf(k) + rsurf if (r2 .le. rbig*rbig) then r = sqrt(r2) rsmall = rpmf(k) - rcarbon part = pisurf * (rbig-r) * (1.0d0+rsmall/r) sasa = sasa * (1.0d0-saterm*part) end if end if end do sasa = acsurf * sasa cutv = tanh(tslope*(sasa-toffset)) cutmtx(i) = 0.5d0 * (1.0d0+cutv) end do ! find the hydrophobic PMF energy via a double loop search do i = 1, n omit(i) = 0 end do do ii = 1, npmf-1 i = ipmf(ii) sschk = 0 do j = 1, n12(i) k = i12(j,i) omit(k) = i if (atomic(k) .eq. 16) sschk = k end do do j = 1, n13(i) k = i13(j,i) omit(k) = i end do do j = 1, n14(i) k = i14(j,i) omit(k) = i if (sschk .ne. 0) then do jj = 1, n12(k) m = i12(jj,k) if (atomic(m) .eq. 16) then do kk = 1, n12(m) if (i12(kk,m) .eq. sschk) omit(k) = 0 end do end if end do end if end do do kk = ii+1, npmf k = ipmf(kk) if (omit(k) .ne. i) then xr = x(i) - x(k) yr = y(i) - y(k) zr = z(i) - z(k) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. hpmfcut2) then r = sqrt(r2) arg1 = (r-c1) * w1 arg12 = arg1 * arg1 arg2 = (r-c2) * w2 arg22 = arg2 * arg2 arg3 = (r-c3) * w3 arg32 = arg3 * arg3 e1 = h1 * exp(-arg12) e2 = h2 * exp(-arg22) e3 = h3 * exp(-arg32) sum = e1 + e2 + e3 e = sum * cutmtx(i) * cutmtx(k) ehp = ehp + e end if end if end do end do ehpmf = ehp return end function ehpmf !+++++++++++++++++++++++++++++++++++++++++++++++++++ ! !"sort" takes an input list of integers and sorts it !into ascending order using the Heapsort algorithm ! !+++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine sort (n,list) implicit none integer i,j,k,n integer index,lists integer list(*) ! perform the heapsort of the input list k = n/2 + 1 index = n do while (n .gt. 1) if (k .gt. 1) then k = k - 1 lists = list(k) else lists = list(index) list(index) = list(1) index = index - 1 if (index .le. 1) then list(1) = lists return end if end if i = k j = k + k do while (j .le. index) if (j .lt. index) then if (list(j) .lt. list(j+1)) j = j + 1 end if if (lists .lt. list(j)) then list(i) = list(j) i = j j = j + j else j = index + 1 end if end do list(i) = lists end do return end subroutine sort end module genhpmf