source/unres/map.f90

   1       module map_
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use names
   5       use geometry_data
   6       use control_data
   7       use energy_data
   8       use map_data
   9       implicit none
  10 !-----------------------------------------------------------------------------
  11 !
  12 !
  13 !-----------------------------------------------------------------------------
  14       contains
  15 !-----------------------------------------------------------------------------
  16 ! check_sc_map.f
  17 !-----------------------------------------------------------------------------
  18       subroutine check_sc_map
  19 ! Subroutine is checking if the fitted function which describs sc_rot_pot
  20 ! is correct, printing, alpha,beta, energy, data - for some known theta.
  21 ! theta angle is read from the input file. Sc_rot_pot are printed
  22 ! for the second  residue in sequance.
  23       use geometry, only:chainbuild
  24       use energy, only:vec_and_deriv,esc
  25 !       include 'DIMENSIONS'
  26 !       include 'COMMON.VAR'
  27 !       include 'COMMON.GEO'
  28 !       include 'COMMON.INTERACT'
  29        real(kind=8) :: xx,yy,zz,al,om
  30        real(kind=8) :: escloc, escloc_min
  31        real(kind=8),dimension(50000) :: escloc_ene,alph_plot,beta_plot
  32        integer,dimension(5000) :: al_plot,be_plot
  33        integer :: iialph, iibet,it
  34
  35 !el local variables
  36        integer :: i,j
  37
  38        write (2,*) "Side-chain-rotamer potential energy map!!!!"
  39        escloc_min = 1000000.00
  40 !       it=itype(2)
  41        i = 0
  42        do iialph=0,18
  43          do iibet=-18,18
  44            i = i + 1
  45            al = iialph*10.0d0*deg2rad
  46            om = iibet*10.0d0*deg2rad
  47            zz = dcos(al)
  48            xx = -dsin(al)*dcos(om)
  49            yy = -dsin(al)*dsin(om)
  50            alph(2)=dacos(xx)
  51            omeg(2)=-datan2(zz,yy)
  52            al_plot(i)=alph(2)*rad2deg
  53            be_plot(i)=omeg(2)*rad2deg
  54 !         write(2,*) alph(2)*rad2deg, omeg(2)*rad2deg
  55            alph_plot(i) = al*rad2deg
  56            beta_plot(i) = om*rad2deg
  57            call chainbuild
  58            call vec_and_deriv
  59            call esc(escloc)
  60            escloc_ene(i) = escloc
  61            if (escloc_min.gt.escloc_ene(i)) escloc_min=escloc_ene(i)
  62          enddo
  63        enddo
  64 !       write (2,*) "escloc_min = ", escloc_min
  65        print *,"i",i
  66        do j = 1,i
  67           write (2,'(3f10.3,2i9,f12.5)') alph_plot(j), &
  68                beta_plot(j),theta(3)*rad2deg, al_plot(j),be_plot(j),&
  69                escloc_ene(j) !- escloc_min
  70        enddo
  71       return
  72       end subroutine check_sc_map
  73 !-----------------------------------------------------------------------------
  74 ! map.f
  75 !-----------------------------------------------------------------------------
  76       subroutine map
  77
  78       use geometry, only:chainbuild,geom_to_var
  79       use energy
  80       use minimm, only:minimize
  81 !      implicit real*8 (a-h,o-z)
  82 !      include 'DIMENSIONS'
  83 !      include 'COMMON.MAP'
  84 !      include 'COMMON.VAR'
  85 !      include 'COMMON.GEO'
  86 !      include 'COMMON.DERIV'
  87 !      include 'COMMON.IOUNITS'
  88 !      include 'COMMON.NAMES'
  89 !      include 'COMMON.CONTROL'
  90 !      include 'COMMON.TORCNSTR'
  91       real(kind=8) :: energia(0:n_ene)
  92       character(len=5) :: angid(4)=reshape((/'PHI  ','THETA','ALPHA',&
  93                                              'OMEGA'/),shape(angid))
  94       real(kind=8) :: ang_list(10)
  95 !el      real(kind=8),dimension(6*nres) :: g,x  !(maxvar)  (maxvar=6*maxres)
  96       real(kind=8),dimension(:),allocatable :: g,x      !(maxvar)  (maxvar=6*maxres)
  97       integer :: nn(10)
  98 !el local variables
  99       integer :: i,iii,ii,j,k,nmax,ntot,nf,iretcode,nfun
 100       real(kind=8) :: etot,gnorm!,fdum
 101       integer,dimension(1) :: uiparm
 102       real(kind=8),dimension(1) :: urparm
 103       allocate(x(6*nres),g(6*nres))
 104
 105       write (iout,'(a,i3,a)')'Energy map constructed in the following ',&
 106              nmap,' groups of variables:'
 107       do i=1,nmap
 108         write (iout,'(2a,i3,a,i3)') angid(kang(i)),' of residues ',&
 109          res1(i),' to ',res2(i)
 110       enddo
 111       nmax=nstep(1)
 112       do i=2,nmap
 113         if (nmax.lt.nstep(i)) nmax=nstep(i)
 114       enddo
 115       ntot=nmax**nmap
 116       iii=0
 117       write (istat,'(1h#,a14,29a15)') (" ",k=1,nmap),&
 118           (ename(print_order(k)),k=1,nprint_ene),"ETOT","GNORM"
 119       do i=0,ntot-1
 120         ii=i
 121         do j=1,nmap
 122           nn(j)=mod(ii,nmax)+1
 123           ii=ii/nmax
 124         enddo
 125         do j=1,nmap
 126           if (nn(j).gt.nstep(j)) goto 10
 127         enddo
 128         iii=iii+1
 129 !d      write (iout,*) i,iii,(nn(j),j=1,nmap)
 130         do j=1,nmap
 131           ang_list(j)=ang_from(j) &
 132              +(nn(j)-1)*(ang_to(j)-ang_from(j))/nstep(j)
 133           do k=res1(j),res2(j)
 134             goto (1,2,3,4), kang(j)
 135     1       phi(k)=deg2rad*ang_list(j)
 136             if (minim) phi0(k-res1(j)+1)=deg2rad*ang_list(j)
 137             goto 5
 138     2       theta(k)=deg2rad*ang_list(j)
 139             goto 5
 140     3       alph(k)=deg2rad*ang_list(j)
 141             goto 5
 142     4       omeg(k)=deg2rad*ang_list(j)
 143     5       continue
 144           enddo ! k
 145         enddo ! j
 146         call chainbuild
 147         if (minim) then
 148          call geom_to_var(nvar,x)
 149          call minimize(etot,x,iretcode,nfun)
 150          print *,'SUMSL return code is',iretcode,' eval ',nfun
 151 !         call intout
 152         else
 153          call zerograd
 154          call geom_to_var(nvar,x)
 155         endif
 156          call etotal(energia)
 157          etot = energia(0)
 158          nf=1
 159          nfl=3
 160          call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
 161          gnorm=0.0d0
 162          do k=1,nvar
 163            gnorm=gnorm+g(k)**2
 164          enddo
 165         etot=energia(0)
 166
 167         gnorm=dsqrt(gnorm)
 168 !        write (iout,'(6(1pe15.5))') (ang_list(k),k=1,nmap),etot,gnorm
 169         write (istat,'(30e15.5)') (ang_list(k),k=1,nmap),&
 170          (energia(print_order(ii)),ii=1,nprint_ene),etot,gnorm
 171 !        write (iout,*) 'POINT',I,' ANGLES:',(ang_list(k),k=1,nmap)
 172 !        call intout
 173 !        call enerprint(energia)
 174    10   continue
 175       enddo ! i
 176       deallocate(x,g)
 177       return
 178       end subroutine map
 179 !-----------------------------------------------------------------------------
 180       subroutine alloc_map_arrays
 181
 182 ! commom.map
 183 !      common /mapp/
 184       allocate(kang(nmap),res1(nmap),res2(nmap),nstep(nmap)) !(maxvar)
 185       allocate(ang_from(nmap),ang_to(nmap)) !(maxvar)
 186
 187       return
 188       end subroutine alloc_map_arrays
 189 !-----------------------------------------------------------------------------
 190 !-----------------------------------------------------------------------------
 191       end module map_