module map_ !----------------------------------------------------------------------------- use io_units use names use geometry_data use control_data use energy_data use map_data implicit none !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !----------------------------------------------------------------------------- ! check_sc_map.f !----------------------------------------------------------------------------- subroutine check_sc_map ! Subroutine is checking if the fitted function which describs sc_rot_pot ! is correct, printing, alpha,beta, energy, data - for some known theta. ! theta angle is read from the input file. Sc_rot_pot are printed ! for the second residue in sequance. use geometry, only:chainbuild use energy, only:vec_and_deriv,esc ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.GEO' ! include 'COMMON.INTERACT' real(kind=8) :: xx,yy,zz,al,om real(kind=8) :: escloc, escloc_min real(kind=8),dimension(50000) :: escloc_ene,alph_plot,beta_plot integer,dimension(5000) :: al_plot,be_plot integer :: iialph, iibet,it !el local variables integer :: i,j write (2,*) "Side-chain-rotamer potential energy map!!!!" escloc_min = 1000000.00 ! it=itype(2) i = 0 do iialph=0,18 do iibet=-18,18 i = i + 1 al = iialph*10.0d0*deg2rad om = iibet*10.0d0*deg2rad zz = dcos(al) xx = -dsin(al)*dcos(om) yy = -dsin(al)*dsin(om) alph(2)=dacos(xx) omeg(2)=-datan2(zz,yy) al_plot(i)=alph(2)*rad2deg be_plot(i)=omeg(2)*rad2deg ! write(2,*) alph(2)*rad2deg, omeg(2)*rad2deg alph_plot(i) = al*rad2deg beta_plot(i) = om*rad2deg call chainbuild call vec_and_deriv call esc(escloc) escloc_ene(i) = escloc if (escloc_min.gt.escloc_ene(i)) escloc_min=escloc_ene(i) enddo enddo ! write (2,*) "escloc_min = ", escloc_min print *,"i",i do j = 1,i write (2,'(3f10.3,2i9,f12.5)') alph_plot(j), & beta_plot(j),theta(3)*rad2deg, al_plot(j),be_plot(j),& escloc_ene(j) !- escloc_min enddo return end subroutine check_sc_map !----------------------------------------------------------------------------- ! map.f !----------------------------------------------------------------------------- subroutine map use geometry, only:chainbuild,geom_to_var use energy use minimm, only:minimize ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.MAP' ! include 'COMMON.VAR' ! include 'COMMON.GEO' ! include 'COMMON.DERIV' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.CONTROL' ! include 'COMMON.TORCNSTR' real(kind=8) :: energia(0:n_ene) character(len=5) :: angid(4)=reshape((/'PHI ','THETA','ALPHA',& 'OMEGA'/),shape(angid)) real(kind=8) :: ang_list(10) !el real(kind=8),dimension(6*nres) :: g,x !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(:),allocatable :: g,x !(maxvar) (maxvar=6*maxres) integer :: nn(10) !el local variables integer :: i,iii,ii,j,k,nmax,ntot,nf,iretcode,nfun real(kind=8) :: etot,gnorm!,fdum integer,dimension(1) :: uiparm real(kind=8),dimension(1) :: urparm allocate(x(6*nres),g(6*nres)) write (iout,'(a,i3,a)')'Energy map constructed in the following ',& nmap,' groups of variables:' do i=1,nmap write (iout,'(2a,i3,a,i3)') angid(kang(i)),' of residues ',& res1(i),' to ',res2(i) enddo nmax=nstep(1) do i=2,nmap if (nmax.lt.nstep(i)) nmax=nstep(i) enddo ntot=nmax**nmap iii=0 write (istat,'(1h#,a14,29a15)') (" ",k=1,nmap),& (ename(print_order(k)),k=1,nprint_ene),"ETOT","GNORM" do i=0,ntot-1 ii=i do j=1,nmap nn(j)=mod(ii,nmax)+1 ii=ii/nmax enddo do j=1,nmap if (nn(j).gt.nstep(j)) goto 10 enddo iii=iii+1 !d write (iout,*) i,iii,(nn(j),j=1,nmap) do j=1,nmap ang_list(j)=ang_from(j) & +(nn(j)-1)*(ang_to(j)-ang_from(j))/nstep(j) do k=res1(j),res2(j) goto (1,2,3,4), kang(j) 1 phi(k)=deg2rad*ang_list(j) if (minim) phi0(k-res1(j)+1)=deg2rad*ang_list(j) goto 5 2 theta(k)=deg2rad*ang_list(j) goto 5 3 alph(k)=deg2rad*ang_list(j) goto 5 4 omeg(k)=deg2rad*ang_list(j) 5 continue enddo ! k enddo ! j call chainbuild if (minim) then call geom_to_var(nvar,x) call minimize(etot,x,iretcode,nfun) print *,'SUMSL return code is',iretcode,' eval ',nfun ! call intout else call zerograd call geom_to_var(nvar,x) endif call etotal(energia) etot = energia(0) nf=1 nfl=3 call gradient(nvar,x,nf,g,uiparm,urparm,fdum) gnorm=0.0d0 do k=1,nvar gnorm=gnorm+g(k)**2 enddo etot=energia(0) gnorm=dsqrt(gnorm) ! write (iout,'(6(1pe15.5))') (ang_list(k),k=1,nmap),etot,gnorm write (istat,'(30e15.5)') (ang_list(k),k=1,nmap),& (energia(print_order(ii)),ii=1,nprint_ene),etot,gnorm ! write (iout,*) 'POINT',I,' ANGLES:',(ang_list(k),k=1,nmap) ! call intout ! call enerprint(energia) 10 continue enddo ! i ! deallocate(x,g) return end subroutine map !----------------------------------------------------------------------------- subroutine alloc_map_arrays ! commom.map ! common /mapp/ allocate(kang(nmap),res1(nmap),res2(nmap),nstep(nmap)) !(maxvar) allocate(ang_from(nmap),ang_to(nmap)) !(maxvar) return end subroutine alloc_map_arrays !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module map_