source/wham/src-M-Pawel/initialize_p.F.org

   1       subroutine initialize
   2 C
   3 C Define constants and zero out tables.
   4 C
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7       include 'DIMENSIONS.ZSCOPT'
   8 #ifdef MPI
   9       include 'mpif.h'
  10 #endif
  11       include 'COMMON.IOUNITS'
  12       include 'COMMON.CHAIN'
  13       include 'COMMON.INTERACT'
  14       include 'COMMON.GEO'
  15       include 'COMMON.LOCAL'
  16       include 'COMMON.TORSION'
  17       include 'COMMON.FFIELD'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.MINIM'
  20       include 'COMMON.DERIV'
  21       include "COMMON.WEIGHTS"
  22       include "COMMON.NAMES"
  23       include "COMMON.TIME1"
  24 C
  25 C The following is just to define auxiliary variables used in angle conversion
  26 C
  27       pi=4.0D0*datan(1.0D0)
  28       dwapi=2.0D0*pi
  29       dwapi3=dwapi/3.0D0
  30       pipol=0.5D0*pi
  31       deg2rad=pi/180.0D0
  32       rad2deg=1.0D0/deg2rad
  33       angmin=10.0D0*deg2rad
  34 C
  35 C Define I/O units.
  36 C
  37       inp=    1
  38       iout=   2
  39       ipdbin= 3
  40       ipdb=   7
  41       imol2=  4
  42       igeom=  8
  43       intin=  9
  44       ithep= 11
  45       irotam=12
  46       itorp= 13
  47       itordp= 23
  48       ielep= 14
  49       isidep=15
  50       iscpp=25
  51       icbase=16
  52       ifourier=20
  53       istat= 17
  54       ientin=18
  55       ientout=19
  56 C
  57 C CSA I/O units (separated from others especially for Jooyoung)
  58 C
  59       icsa_rbank=30
  60       icsa_seed=31
  61       icsa_history=32
  62       icsa_bank=33
  63       icsa_bank1=34
  64       icsa_alpha=35
  65       icsa_alpha1=36
  66       icsa_bankt=37
  67       icsa_int=39
  68       icsa_bank_reminimized=38
  69       icsa_native_int=41
  70       icsa_in=40
  71 C
  72 C Set default weights of the energy terms.
  73 C
  74       wlong=1.0D0
  75       welec=1.0D0
  76       wtor =1.0D0
  77       wang =1.0D0
  78       wscloc=1.0D0
  79       wstrain=1.0D0
  80 C
  81 C Zero out tables.
  82 C
  83       ndih_constr=0
  84       do i=1,maxres2
  85         do j=1,3
  86           c(j,i)=0.0D0
  87           dc(j,i)=0.0D0
  88         enddo
  89       enddo
  90       do i=1,maxres
  91         do j=1,3
  92           xloc(j,i)=0.0D0
  93         enddo
  94       enddo
  95       do i=1,ntyp
  96         do j=1,ntyp
  97           aa(i,j)=0.0D0
  98           bb(i,j)=0.0D0
  99           augm(i,j)=0.0D0
 100           sigma(i,j)=0.0D0
 101           r0(i,j)=0.0D0
 102           chi(i,j)=0.0D0
 103         enddo
 104         do j=1,2
 105           bad(i,j)=0.0D0
 106         enddo
 107         chip(i)=0.0D0
 108         alp(i)=0.0D0
 109         sigma0(i)=0.0D0
 110         sigii(i)=0.0D0
 111         rr0(i)=0.0D0
 112         a0thet(i)=0.0D0
 113         do j=1,2
 114           athet(j,i)=0.0D0
 115           bthet(j,i)=0.0D0
 116         enddo
 117         do j=0,3
 118           polthet(j,i)=0.0D0
 119         enddo
 120         do j=1,3
 121           gthet(j,i)=0.0D0
 122         enddo
 123         theta0(i)=0.0D0
 124         sig0(i)=0.0D0
 125         sigc0(i)=0.0D0
 126         do j=1,maxlob
 127           bsc(j,i)=0.0D0
 128           do k=1,3
 129             censc(k,j,i)=0.0D0
 130           enddo
 131           do k=1,3
 132             do l=1,3
 133               gaussc(l,k,j,i)=0.0D0
 134             enddo
 135           enddo
 136           nlob(i)=0
 137         enddo
 138       enddo
 139       nlob(ntyp1)=0
 140       dsc(ntyp1)=0.0D0
 141       do i=1,maxtor
 142         itortyp(i)=0
 143         do j=1,maxtor
 144           do k=1,maxterm
 145             v1(k,j,i)=0.0D0
 146             v2(k,j,i)=0.0D0
 147           enddo
 148         enddo
 149       enddo
 150       do i=1,maxres
 151         itype(i)=0
 152         itel(i)=0
 153       enddo
 154 C Initialize the bridge arrays
 155       ns=0
 156       nss=0
 157       nhpb=0
 158       do i=1,maxss
 159         iss(i)=0
 160       enddo
 161       do i=1,maxdim
 162         dhpb(i)=0.0D0
 163       enddo
 164       do i=1,maxres
 165         ihpb(i)=0
 166         jhpb(i)=0
 167       enddo
 168 C
 169 C Initialize timing.
 170 C
 171       call set_timers
 172 C
 173 C Initialize variables used in minimization.
 174 C
 175 c     maxfun=5000
 176 c     maxit=2000
 177       maxfun=500
 178       maxit=200
 179       tolf=1.0D-2
 180       rtolf=5.0D-4
 181 C
 182 C Initialize the variables responsible for the mode of gradient storage.
 183 C
 184       nfl=0
 185       icg=1
 186       do i=1,14
 187         do j=1,14
 188           if (print_order(i).eq.j) then
 189             iw(print_order(i))=j
 190             goto 1121
 191           endif
 192         enddo
 193 1121    continue
 194       enddo
 195       calc_grad=.false.
 196 C Set timers and counters for the respective routines
 197       t_func = 0.0d0
 198       t_grad = 0.0d0
 199       t_fhel = 0.0d0
 200       t_fbet = 0.0d0
 201       t_ghel = 0.0d0
 202       t_gbet = 0.0d0
 203       t_viol = 0.0d0
 204       t_gviol = 0.0d0
 205       n_func = 0
 206       n_grad = 0
 207       n_fhel = 0
 208       n_fbet = 0
 209       n_ghel = 0
 210       n_gbet = 0
 211       n_viol = 0
 212       n_gviol = 0
 213       n_map = 0
 214       return
 215       end
 216 c-------------------------------------------------------------------------
 217       block data nazwy
 218       implicit real*8 (a-h,o-z)
 219       include 'DIMENSIONS'
 220       include 'DIMENSIONS.ZSCOPT'
 221       include 'COMMON.NAMES'
 222       include 'COMMON.FFIELD'
 223       data restyp /
 224      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 225      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
 226       data onelet /
 227      &'C','M','F','I','L','V','W','Y','A','G','T',
 228      &'S','Q','N','E','D','H','R','K','P','X'/
 229       data potname /'LJ','LJK','BP','GB','GBV'/
 230       data ename /
 231      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 232      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 233      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EVDW2_14",2*" "/
 234       data wname /
 235      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 236      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 237      &   "SCAL14",2*" "/
 238 #ifdef SCP14
 239       data nprint_ene /15/
 240       data print_order /1,2,3,11,12,13,14,4,5,6,7,8,9,10,16,0/
 241 #else
 242       data nprint_ene /14/
 243       data print_order /1,2,3,11,12,13,14,4,5,6,7,8,9,10,3*0/
 244 #endif
 245       end
 246 c---------------------------------------------------------------------------
 247       subroutine init_int_table
 248       implicit real*8 (a-h,o-z)
 249       include 'DIMENSIONS'
 250       include 'DIMENSIONS.ZSCOPT'
 251 #ifdef MPI
 252       include 'mpif.h'
 253 #endif
 254 #ifdef MP
 255       include 'COMMON.INFO'
 256 #endif
 257       include 'COMMON.CHAIN'
 258       include 'COMMON.INTERACT'
 259       include 'COMMON.LOCAL'
 260       include 'COMMON.SBRIDGE'
 261       include 'COMMON.IOUNITS'
 262       logical scheck,lprint
 263 #ifdef MPL
 264       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 265      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 266 C... Determine the numbers of start and end SC-SC interaction
 267 C... to deal with by current processor.
 268       lprint=.false.
 269       if (lprint)
 270      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 271       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 272       MyRank=MyID-(MyGroup-1)*fgProcs
 273       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 274       if (lprint)
 275      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 276      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 277      &  ' my_sc_inde',my_sc_inde
 278       ind_sctint=0
 279       iatsc_s=0
 280       iatsc_e=0
 281 #endif
 282       lprint=.false.
 283       do i=1,maxres
 284         nint_gr(i)=0
 285         nscp_gr(i)=0
 286         do j=1,maxint_gr
 287           istart(i,1)=0
 288           iend(i,1)=0
 289           ielstart(i)=0
 290           ielend(i)=0
 291           iscpstart(i,1)=0
 292           iscpend(i,1)=0
 293         enddo
 294       enddo
 295       ind_scint=0
 296       ind_scint_old=0
 297 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 298 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 299       do i=nnt,nct-1
 300         scheck=.false.
 301         do ii=1,nss
 302           if (ihpb(ii).eq.i+nres) then
 303             scheck=.true.
 304             jj=jhpb(ii)-nres
 305             goto 10
 306           endif
 307         enddo
 308    10   continue
 309 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 310         if (scheck) then
 311           if (jj.eq.i+1) then
 312 #ifdef MPL
 313             write (iout,*) 'jj=i+1'
 314             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 315      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 316 #else
 317             nint_gr(i)=1
 318             istart(i,1)=i+2
 319             iend(i,1)=nct
 320 #endif
 321           else if (jj.eq.nct) then
 322 #ifdef MPL
 323             write (iout,*) 'jj=nct'
 324             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 325      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 326 #else
 327             nint_gr(i)=1
 328             istart(i,1)=i+1
 329             iend(i,1)=nct-1
 330 #endif
 331           else
 332 #ifdef MPL
 333             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 334      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 335             ii=nint_gr(i)+1
 336             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 337      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 338 #else
 339             nint_gr(i)=2
 340             istart(i,1)=i+1
 341             iend(i,1)=jj-1
 342             istart(i,2)=jj+1
 343             iend(i,2)=nct
 344 #endif
 345           endif
 346         else
 347 #ifdef MPL
 348           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 349      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 350 #else
 351           nint_gr(i)=1
 352           istart(i,1)=i+1
 353           iend(i,1)=nct
 354           ind_scint=int_scint+nct-i
 355 #endif
 356         endif
 357 #ifdef MPL
 358         ind_scint_old=ind_scint
 359 #endif
 360       enddo
 361    12 continue
 362 #ifndef MPL
 363       iatsc_s=nnt
 364       iatsc_e=nct-1
 365 #endif
 366 #ifdef MPL
 367       if (lprint) then
 368         write (iout,*) 'Processor',MyID,' Group',MyGroup
 369         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 370       endif
 371 #endif
 372       if (lprint) then
 373       write (iout,'(a)') 'Interaction array:'
 374       do i=iatsc_s,iatsc_e
 375         write (iout,'(i3,2(2x,2i3))')
 376      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 377       enddo
 378       endif
 379       ispp=2
 380 #ifdef MPL
 381 C Now partition the electrostatic-interaction array
 382       npept=nct-nnt
 383       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 384       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 385       if (lprint)
 386      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 387      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 388      &               ' my_ele_inde',my_ele_inde
 389       iatel_s=0
 390       iatel_e=0
 391       ind_eleint=0
 392       ind_eleint_old=0
 393       do i=nnt,nct-3
 394         ijunk=0
 395         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 396      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 397       enddo ! i
 398    13 continue
 399 #else
 400       iatel_s=nnt
 401       iatel_e=nct-3
 402       do i=iatel_s,iatel_e
 403         ielstart(i)=i+2
 404         ielend(i)=nct-1
 405       enddo
 406 #endif
 407       if (lprint) then
 408         write (iout,'(a)') 'Electrostatic interaction array:'
 409         do i=iatel_s,iatel_e
 410           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 411         enddo
 412       endif ! lprint
 413 c     iscp=3
 414       iscp=2
 415 C Partition the SC-p interaction array
 416 #ifdef MPL
 417       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 418       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 419       if (lprint)
 420      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 421      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 422      &               ' my_scp_inde',my_scp_inde
 423       iatscp_s=0
 424       iatscp_e=0
 425       ind_scpint=0
 426       ind_scpint_old=0
 427       do i=nnt,nct-1
 428         if (i.lt.nnt+iscp) then
 429 cd        write (iout,*) 'i.le.nnt+iscp'
 430           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 431      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 432      &      iscpend(i,1),*14)
 433         else if (i.gt.nct-iscp) then
 434 cd        write (iout,*) 'i.gt.nct-iscp'
 435           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 436      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 437      &      iscpend(i,1),*14)
 438         else
 439           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 440      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 441      &      iscpend(i,1),*14)
 442           ii=nscp_gr(i)+1
 443           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 444      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 445      &      iscpend(i,ii),*14)
 446         endif
 447       enddo ! i
 448    14 continue
 449 #else
 450       iatscp_s=nnt
 451       iatscp_e=nct-1
 452       do i=nnt,nct-1
 453         if (i.lt.nnt+iscp) then
 454           nscp_gr(i)=1
 455           iscpstart(i,1)=i+iscp
 456           iscpend(i,1)=nct
 457         elseif (i.gt.nct-iscp) then
 458           nscp_gr(i)=1
 459           iscpstart(i,1)=nnt
 460           iscpend(i,1)=i-iscp
 461         else
 462           nscp_gr(i)=2
 463           iscpstart(i,1)=nnt
 464           iscpend(i,1)=i-iscp
 465           iscpstart(i,2)=i+iscp
 466           iscpend(i,2)=nct
 467         endif
 468       enddo ! i
 469 #endif
 470       if (lprint) then
 471         write (iout,'(a)') 'SC-p interaction array:'
 472         do i=iatscp_s,iatscp_e
 473           write (iout,'(i3,2(2x,2i3))')
 474      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 475         enddo
 476       endif ! lprint
 477 C Partition local interactions
 478 #ifdef MPL
 479       call int_bounds(nres-2,loc_start,loc_end)
 480       loc_start=loc_start+1
 481       loc_end=loc_end+1
 482       call int_bounds(nres-2,ithet_start,ithet_end)
 483       ithet_start=ithet_start+2
 484       ithet_end=ithet_end+2
 485       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 486       iphi_start=iphi_start+nnt+2
 487       iphi_end=iphi_end+nnt+2
 488       if (lprint) then
 489         write (iout,*) 'Processor:',MyID,
 490      & ' loc_start',loc_start,' loc_end',loc_end,
 491      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 492      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 493         write (*,*) 'Processor:',MyID,
 494      & ' loc_start',loc_start,' loc_end',loc_end,
 495      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 496      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 497       endif
 498       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 499         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 500      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 501      & ' SC-p interactions','were distributed among',fgprocs,
 502      & ' fine-grain processors.'
 503       endif
 504 #else
 505       loc_start=2
 506       loc_end=nres-1
 507       ithet_start=3
 508       ithet_end=nres
 509       iphi_start=nnt+3
 510       iphi_end=nct
 511 #endif
 512       return
 513       end
 514 c---------------------------------------------------------------------------
 515       subroutine int_partition(int_index,lower_index,upper_index,atom,
 516      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
 517       implicit real*8 (a-h,o-z)
 518       include 'DIMENSIONS'
 519       include 'COMMON.IOUNITS'
 520       integer int_index,lower_index,upper_index,atom,at_start,at_end,
 521      & first_atom,last_atom,int_gr,jat_start,jat_end
 522       logical lprn
 523       lprn=.false.
 524       if (lprn) write (iout,*) 'int_index=',int_index
 525       int_index_old=int_index
 526       int_index=int_index+last_atom-first_atom+1
 527       if (lprn)
 528      &   write (iout,*) 'int_index=',int_index,
 529      &               ' int_index_old',int_index_old,
 530      &               ' lower_index=',lower_index,
 531      &               ' upper_index=',upper_index,
 532      &               ' atom=',atom,' first_atom=',first_atom,
 533      &               ' last_atom=',last_atom
 534       if (int_index.ge.lower_index) then
 535         int_gr=int_gr+1
 536         if (at_start.eq.0) then
 537           at_start=atom
 538           jat_start=first_atom-1+lower_index-int_index_old
 539         else
 540           jat_start=first_atom
 541         endif
 542         if (lprn) write (iout,*) 'jat_start',jat_start
 543         if (int_index.ge.upper_index) then
 544           at_end=atom
 545           jat_end=first_atom-1+upper_index-int_index_old
 546           return1
 547         else
 548           jat_end=last_atom
 549         endif
 550         if (lprn) write (iout,*) 'jat_end',jat_end
 551       endif
 552       return
 553       end
 554 c------------------------------------------------------------------------------
 555       subroutine hpb_partition
 556       implicit real*8 (a-h,o-z)
 557       include 'DIMENSIONS'
 558       include 'COMMON.SBRIDGE'
 559       include 'COMMON.IOUNITS'
 560 #ifdef MPL
 561       include 'COMMON.INFO'
 562       call int_bounds(nhpb,link_start,link_end)
 563 #else
 564       link_start=1
 565       link_end=nhpb
 566 #endif
 567 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 568 cd   &  ' nhpb',nhpb,' link_start=',link_start,
 569 cd   &  ' link_end',link_end
 570       return
 571       end