source/cluster/wham/src-M/initialize_p.F

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