source/unres/src-HCD-5D/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit none
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.SPLITELE'
  12       integer i,j
  13       logical file_exist
  14 C Read job setup parameters
  15       call read_control
  16 C Read force-field parameters except weights
  17       call parmread
  18 C Read control parameters for energy minimzation if required
  19       if (minim) call read_minim
  20 C Read MCM control parameters if required
  21       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  22 C Read MD control parameters if reqjuired
  23       if (modecalc.eq.12) call read_MDpar
  24 C Read MREMD control parameters if required
  25       if (modecalc.eq.14) then
  26          call read_MDpar
  27          call read_REMDpar
  28       endif
  29 C Read MUCA control parameters if required
  30       if (lmuca) call read_muca
  31 C Read CSA control parameters if required (from fort.40 if exists
  32 C otherwise from general input file)
  33       if (modecalc.eq.8) then
  34        inquire (file="fort.40",exist=file_exist)
  35        if (.not.file_exist) call csaread
  36       endif
  37 cfmc      if (modecalc.eq.10) call mcmfread
  38 c Read energy-term weights and disulfide parameters
  39       call weightread
  40 C Read molecule information, molecule geometry, energy-term weights, and
  41 C restraints if requested
  42       call molread
  43 C Print restraint information
  44 #ifdef MPI
  45       if (.not. out1file .or. me.eq.king) then
  46 #endif
  47       if (nhpb.gt.nss)  then
  48         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  49      &  "The following",nhpb-nss,
  50      &  " distance restraints have been imposed:",
  51      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  52      &  "  score"," type"
  53         do i=nss+1,nhpb
  54           write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
  55      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
  56      &  irestr_type(i)
  57         enddo
  58       endif
  59       if (npeak.gt.0) then
  60         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  61      &  "The following",npeak,
  62      &  " NMR peak restraints have been imposed:",
  63      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  64      &  "  score"," type"," ipeak"
  65         do i=1,npeak
  66           do j=ipeak(1,i),ipeak(2,i)
  67             write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
  68      &       jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
  69      &       forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
  70           enddo
  71         enddo
  72         write (iout,*) "The ipeak array"
  73         do i=1,npeak
  74           write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
  75         enddo
  76       endif
  77 #ifdef MPI
  78       endif
  79 #endif
  80 c      print *,"Processor",myrank," leaves READRTNS"
  81       return
  82       end
  83 C-------------------------------------------------------------------------------
  84       subroutine read_control
  85 C
  86 C Read contorl data
  87 C
  88       implicit none
  89       include 'DIMENSIONS'
  90 #ifdef MP
  91       include 'mpif.h'
  92       logical OKRandom, prng_restart
  93       double precision  r1
  94 #endif
  95       include 'COMMON.IOUNITS'
  96       include 'COMMON.TIME1'
  97       include 'COMMON.THREAD'
  98       include 'COMMON.SBRIDGE'
  99       include 'COMMON.CONTROL'
 100       include 'COMMON.SAXS'
 101       include 'COMMON.MCM'
 102       include 'COMMON.MAP'
 103       include 'COMMON.HEADER'
 104       include 'COMMON.CSA'
 105       include 'COMMON.CHAIN'
 106       include 'COMMON.MUCA'
 107       include 'COMMON.MD'
 108       include 'COMMON.FFIELD'
 109       include 'COMMON.INTERACT'
 110       include 'COMMON.SETUP'
 111       include 'COMMON.SPLITELE'
 112       include 'COMMON.SHIELD'
 113       include 'COMMON.GEO'
 114       integer i
 115       integer KDIAG,ICORFL,IXDR
 116       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
 117       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
 118       character*80 ucase
 119       character*320 controlcard
 120       double precision seed
 121
 122       nglob_csa=0
 123       eglob_csa=1d99
 124       nmin_csa=0
 125       read (INP,'(a)') titel
 126       call card_concat(controlcard)
 127 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
 128 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 129       call reada(controlcard,'SEED',seed,0.0D0)
 130       call random_init(seed)
 131 C Set up the time limit (caution! The time must be input in minutes!)
 132       read_cart=index(controlcard,'READ_CART').gt.0
 133       out_cart=index(controlcard,'OUT_CART').gt.0
 134       out_int=index(controlcard,'OUT_INT').gt.0
 135       gmatout=index(controlcard,'GMATOUT').gt.0
 136       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 137 C this variable with_theta_constr is the variable which allow to read and execute the
 138 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 139       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 140       write (iout,*) "with_theta_constr ",with_theta_constr
 141       call readi(controlcard,'NSAXS',nsaxs,0)
 142       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 143       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 144       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 145       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 146      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 147       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
 148       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
 149       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
 150       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
 151       call readi(controlcard,'SYM',symetr,1)
 152       call readi(controlcard,'PERMUT',npermut,1)
 153       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
 154       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 155       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 156 c      call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 157 c      call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 158 c      call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 159 c      call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 160 c      call reada(controlcard,'DRMS',drms,0.1D0)
 161 c      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 162 c       write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 163 c       write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 164 c       write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 165 c       write (iout,'(a,f10.1)')'DRMS    = ',drms
 166 cc       write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 167 c       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 168 c      endif
 169       call readi(controlcard,'NZ_START',nz_start,0)
 170       call readi(controlcard,'NZ_END',nz_end,0)
 171       call readi(controlcard,'NRAN_START',nran_start,0)
 172       write (iout,*) "nran_start",nran_start
 173 c      call readi(controlcard,'IZ_SC',iz_sc,0)
 174       timlim=60.0D0*timlim
 175       safety = 60.0d0*safety
 176       modecalc=0
 177       call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
 178       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 179       minim=(index(controlcard,'MINIMIZE').gt.0)
 180       dccart=(index(controlcard,'CART').gt.0)
 181       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 182       overlapsc=.not.overlapsc
 183       searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
 184      &  index(controlcard,'NOSEARCHSC').eq.0)
 185       sideadd=(index(controlcard,'SIDEADD').gt.0)
 186       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 187       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
 188       outpdb=(index(controlcard,'PDBOUT').gt.0)
 189       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 190       pdbref=(index(controlcard,'PDBREF').gt.0)
 191       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 192       indpdb=index(controlcard,'PDBSTART')
 193       AFMlog=(index(controlcard,'AFM'))
 194       selfguide=(index(controlcard,'SELFGUIDE'))
 195 c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
 196       call readi(controlcard,'TUBEMOD',tubelog,0)
 197 c      write (iout,*) TUBElog,"TUBEMODE"
 198       call readi(controlcard,'IPRINT',iprint,0)
 199 c 6/22/2021 AL: alpha_GB: parameter to switch between the GB SC-SC
 200 c    interaction potential and the all-repulsive potential with singularity
 201 c    at zero site-site distance
 202       call reada(controlcard,'ALPHA_GB',alpha_GB,1.0d-2)
 203       alpha_GB1 = 1.0d0+1.0d0/alpha_GB
 204       write (iout,*) "alpha_GB",alpha_GB," alpha_GB1",alpha_GB1
 205 C SHIELD keyword sets if the shielding effect of side-chains is used
 206 C 0 denots no shielding is used all peptide are equally despite the
 207 C solvent accesible area
 208 C 1 the newly introduced function
 209 C 2 reseved for further possible developement
 210       call readi(controlcard,'SHIELD',shield_mode,0)
 211 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 212         write(iout,*) "shield_mode",shield_mode
 213 C      endif
 214       call readi(controlcard,'TORMODE',tor_mode,0)
 215 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 216         write(iout,*) "torsional and valence angle mode",tor_mode
 217       call readi(controlcard,'MAXGEN',maxgen,10000)
 218       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 219       call readi(controlcard,"KDIAG",kdiag,0)
 220       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 221       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 222      & write (iout,*) "RESCALE_MODE",rescale_mode
 223       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 224       if (index(controlcard,'REGULAR').gt.0.0D0) then
 225         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 226         modecalc=1
 227         refstr=.true.
 228       endif
 229       if (refstr) then
 230       if (index(controlcard,"CASC").gt.0) then
 231         iz_sc=1
 232 c      else if (index(controlcard,"CAONLY").gt.0) then
 233 c        iz_sc=0
 234       else if (index(controlcard,"SCONLY").gt.0) then
 235         iz_sc=2
 236       else
 237         iz_sc=0
 238 c        write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
 239 c     &    " specifying REFSTR, PDBREF or REGULAR."
 240 c        stop
 241       endif
 242       endif
 243       if (index(controlcard,'CHECKGRAD').gt.0) then
 244         modecalc=5
 245         if (index(controlcard,' CART').gt.0) then
 246           icheckgrad=1
 247         elseif (index(controlcard,'CARINT').gt.0) then
 248           icheckgrad=2
 249         else
 250           icheckgrad=3
 251         endif
 252         call reada(controlcard,'DELTA',aincr,1.0d-7)
 253 c        write (iout,*) "icheckgrad",icheckgrad
 254       elseif (index(controlcard,'THREAD').gt.0) then
 255         modecalc=2
 256         call readi(controlcard,'THREAD',nthread,0)
 257         if (nthread.gt.0) then
 258           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 259         else
 260           if (fg_rank.eq.0)
 261      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 262           stop 'Error termination in Read_Control.'
 263         endif
 264       else if (index(controlcard,'MCMA').gt.0) then
 265         modecalc=3
 266       else if (index(controlcard,'MCEE').gt.0) then
 267         modecalc=6
 268       else if (index(controlcard,'MULTCONF').gt.0) then
 269         modecalc=4
 270       else if (index(controlcard,'MAP').gt.0) then
 271         modecalc=7
 272         call readi(controlcard,'MAP',nmap,0)
 273       else if (index(controlcard,'CSA').gt.0) then
 274         modecalc=8
 275 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 276 crc
 277 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 278 crc
 279 crc        modecalc=9
 280 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 281 cfmc        modecalc=10
 282       else if (index(controlcard,'SOFTREG').gt.0) then
 283         modecalc=11
 284       else if (index(controlcard,'CHECK_BOND').gt.0) then
 285         modecalc=-1
 286       else if (index(controlcard,'TEST').gt.0) then
 287         modecalc=-2
 288       else if (index(controlcard,'MD').gt.0) then
 289         modecalc=12
 290       else if (index(controlcard,'RE ').gt.0) then
 291         modecalc=14
 292       endif
 293
 294       lmuca=index(controlcard,'MUCA').gt.0
 295       call readi(controlcard,'MUCADYN',mucadyn,0)
 296       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 297       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 298      & then
 299        write (iout,*) 'MUCADYN=',mucadyn
 300        write (iout,*) 'MUCASMOOTH=',muca_smooth
 301       endif
 302
 303       iscode=index(controlcard,'ONE_LETTER')
 304       indphi=index(controlcard,'PHI')
 305       indback=index(controlcard,'BACK')
 306       iranconf=index(controlcard,'RAND_CONF')
 307       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
 308       extconf=(index(controlcard,'EXTCONF').gt.0)
 309       if (start_from_model) then
 310         iranconf=0
 311         extconf=.false.
 312       endif
 313       i2ndstr=index(controlcard,'USE_SEC_PRED')
 314       gradout=index(controlcard,'GRADOUT').gt.0
 315       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 316 C  DISTCHAINMAX become obsolete for periodic boundry condition
 317       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 318 C Reading the dimensions of box in x,y,z coordinates
 319       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 320       call reada(controlcard,'BOXY',boxysize,100.0d0)
 321       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 322       write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
 323 c Cutoff range for interactions
 324       call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
 325       call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
 326       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 327       write (iout,*) "Cutoff on interactions",r_cut_int
 328       write (iout,*)
 329      & "Cutoff in switching short and long range interactions in RESPA",
 330      & r_cut_respa
 331       write (iout,*) "lambda in switch function",rlamb
 332       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 333       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 334       if (lipthick.gt.0.0d0) then
 335        bordliptop=(boxzsize+lipthick)/2.0
 336        bordlipbot=bordliptop-lipthick
 337 C      endif
 338       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 339      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 340       buflipbot=bordlipbot+lipbufthick
 341       bufliptop=bordliptop-lipbufthick
 342       if ((lipbufthick*2.0d0).gt.lipthick)
 343      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 344       endif
 345       write(iout,*) "bordliptop=",bordliptop
 346       write(iout,*) "bordlipbot=",bordlipbot
 347       write(iout,*) "bufliptop=",bufliptop
 348       write(iout,*) "buflipbot=",buflipbot
 349       write (iout,*) "SHIELD MODE",shield_mode
 350       if (TUBElog.gt.0) then
 351        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 352        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 353        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 354        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 355        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 356        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 357        buftubebot=bordtubebot+tubebufthick
 358        buftubetop=bordtubetop-tubebufthick
 359       endif
 360       if (me.eq.king .or. .not.out1file )
 361      &  write (iout,*) "DISTCHAINMAX",distchainmax
 362
 363       if(me.eq.king.or..not.out1file)
 364      & write (iout,'(2a)') diagmeth(kdiag),
 365      &  ' routine used to diagonalize matrices.'
 366       return
 367       end
 368 c--------------------------------------------------------------------------
 369       subroutine read_REMDpar
 370 C
 371 C Read REMD settings
 372 C
 373       implicit none
 374       include 'DIMENSIONS'
 375       include 'COMMON.IOUNITS'
 376       include 'COMMON.TIME1'
 377       include 'COMMON.MD'
 378 #ifndef LANG0
 379       include 'COMMON.LANGEVIN'
 380 #else
 381 #ifdef FIVEDIAG
 382       include 'COMMON.LANGEVIN.lang0.5diag'
 383 #else
 384       include 'COMMON.LANGEVIN.lang0'
 385 #endif
 386 #endif
 387       include 'COMMON.INTERACT'
 388       include 'COMMON.NAMES'
 389       include 'COMMON.GEO'
 390       include 'COMMON.REMD'
 391       include 'COMMON.CONTROL'
 392       include 'COMMON.SETUP'
 393       character*80 ucase
 394       character*320 controlcard
 395       character*3200 controlcard1
 396       integer iremd_m_total,i
 397
 398       if(me.eq.king.or..not.out1file)
 399      & write (iout,*) "REMD setup"
 400
 401       call card_concat(controlcard)
 402       call readi(controlcard,"NREP",nrep,3)
 403       call readi(controlcard,"NSTEX",nstex,1000)
 404       call reada(controlcard,"RETMIN",retmin,10.0d0)
 405       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 406       mremdsync=(index(controlcard,'SYNC').gt.0)
 407       call readi(controlcard,"NSYN",i_sync_step,100)
 408       restart1file=(index(controlcard,'REST1FILE').gt.0)
 409       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 410       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 411       if(max_cache_traj_use.gt.max_cache_traj)
 412      &           max_cache_traj_use=max_cache_traj
 413       if(me.eq.king.or..not.out1file) then
 414 cd       if (traj1file) then
 415 crc caching is in testing - NTWX is not ignored
 416 cd        write (iout,*) "NTWX value is ignored"
 417 cd        write (iout,*) "  trajectory is stored to one file by master"
 418 cd        write (iout,*) "  before exchange at NSTEX intervals"
 419 cd       endif
 420        write (iout,*) "NREP= ",nrep
 421        write (iout,*) "NSTEX= ",nstex
 422        write (iout,*) "SYNC= ",mremdsync
 423        write (iout,*) "NSYN= ",i_sync_step
 424        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 425       endif
 426       remd_tlist=.false.
 427       if (index(controlcard,'TLIST').gt.0) then
 428          remd_tlist=.true.
 429          call card_concat(controlcard1)
 430          read(controlcard1,*) (remd_t(i),i=1,nrep)
 431          if(me.eq.king.or..not.out1file)
 432      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 433       endif
 434       remd_mlist=.false.
 435       if (index(controlcard,'MLIST').gt.0) then
 436          remd_mlist=.true.
 437          call card_concat(controlcard1)
 438          read(controlcard1,*) (remd_m(i),i=1,nrep)
 439          if(me.eq.king.or..not.out1file) then
 440           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 441           iremd_m_total=0
 442           do i=1,nrep
 443            iremd_m_total=iremd_m_total+remd_m(i)
 444           enddo
 445            write (iout,*) 'Total number of replicas ',iremd_m_total
 446           endif
 447          endif
 448         if (adaptive) then
 449           write (iout,*)
 450      &      "Adaptive (PMF-biased) umbrella sampling will be run"
 451           call PMFread
 452         endif
 453       if(me.eq.king.or..not.out1file)
 454      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 455       return
 456       end
 457 c--------------------------------------------------------------------------
 458       subroutine read_MDpar
 459 C
 460 C Read MD settings
 461 C
 462       implicit none
 463       include 'DIMENSIONS'
 464       include 'COMMON.IOUNITS'
 465       include 'COMMON.TIME1'
 466       include 'COMMON.MD'
 467       include 'COMMON.QRESTR'
 468 #ifndef LANG0
 469       include 'COMMON.LANGEVIN'
 470 #else
 471 #ifdef FIVEDIAG
 472       include 'COMMON.LANGEVIN.lang0.5diag'
 473 #else
 474       include 'COMMON.LANGEVIN.lang0'
 475 #endif
 476 #endif
 477       include 'COMMON.INTERACT'
 478       include 'COMMON.NAMES'
 479       include 'COMMON.GEO'
 480       include 'COMMON.SETUP'
 481       include 'COMMON.CONTROL'
 482       include 'COMMON.SPLITELE'
 483       include 'COMMON.FFIELD'
 484       character*80 ucase
 485       character*320 controlcard
 486       integer i
 487       double precision eta
 488
 489       call card_concat(controlcard)
 490       call readi(controlcard,"NSTEP",n_timestep,1000000)
 491       call readi(controlcard,"NTWE",ntwe,100)
 492       call readi(controlcard,"NTWX",ntwx,1000)
 493       call readi(controlcard,"REST_FREQ",irest_freq,1000)
 494       call reada(controlcard,"DT",d_time,1.0d-1)
 495       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 496       call reada(controlcard,"DAMAX",damax,1.0d1)
 497       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 498       call readi(controlcard,"LANG",lang,0)
 499       RESPA = index(controlcard,"RESPA") .gt. 0
 500       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 501       ntime_split0=ntime_split
 502       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 503       ntime_split0=ntime_split
 504 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 505 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 506       rest = index(controlcard,"REST").gt.0
 507       tbf = index(controlcard,"TBF").gt.0
 508       usampl = index(controlcard,"USAMPL").gt.0
 509       scale_umb =  index(controlcard,"SCALE_UMB").gt.0
 510       adaptive = index(controlcard,"ADAPTIVE").gt.0
 511       mdpdb = index(controlcard,"MDPDB").gt.0
 512       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 513       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 514       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 515       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 516       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 517       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 518       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 519       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 520       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 521       large = index(controlcard,"LARGE").gt.0
 522       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 523       rattle = index(controlcard,"RATTLE").gt.0
 524       preminim = index(controlcard,"PREMINIM").gt.0
 525       if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
 526         preminim=.true.
 527         dccart=.true.
 528       endif
 529       write (iout,*) "PREMINIM ",preminim
 530       if (preminim) then
 531         if (index(controlcard,'CART').gt.0) dccart=.true.
 532         write (iout,*) "dccart ",dccart
 533         write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
 534         if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
 535       endif
 536 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 537       nset=0
 538       if(usampl) then
 539         write (iout,*) "Umbrella sampling will be run"
 540         if (scale_umb.and.adaptive) then
 541           write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
 542           write (iout,*) "Select one of those and re-run the job."
 543           stop
 544         endif
 545         if (scale_umb) write (iout,*)
 546      &"Umbrella-restraint force constants will be scaled by temperature"
 547         call read_fragments
 548       endif
 549
 550       if(me.eq.king.or..not.out1file) then
 551        write (iout,*)
 552        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 553        write (iout,*)
 554        write (iout,'(a)') "The units are:"
 555        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 556        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 557      &  " acceleration: angstrom/(48.9 fs)**2"
 558        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 559        write (iout,*)
 560        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 561        write (iout,'(a60,f10.5,a)')
 562      &  "Initial time step of numerical integration:",d_time,
 563      &  " natural units"
 564        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 565        write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
 566      &   " A"
 567        write(iout,'(a60,i5)')"Frequency of updating interaction list",
 568      &   imatupdate
 569        write(iout,'(a60,i5)')"Restart writing frequency",irest_freq
 570        if (RESPA) then
 571         write (iout,'(2a,i4,a)')
 572      &    "A-MTS algorithm used; initial time step for fast-varying",
 573      &    " short-range forces split into",ntime_split," steps."
 574         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 575      &   r_cut_respa," lambda",rlamb
 576        endif
 577        write (iout,'(2a,f10.5)')
 578      &  "Maximum acceleration threshold to reduce the time step",
 579      &  "/increase split number:",damax
 580        write (iout,'(2a,f10.5)')
 581      &  "Maximum predicted energy drift to reduce the timestep",
 582      &  "/increase split number:",edriftmax
 583        write (iout,'(a60,f10.5)')
 584      & "Maximum velocity threshold to reduce velocities:",dvmax
 585        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 586        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 587        if (rattle) write (iout,'(a60)')
 588      &  "Rattle algorithm used to constrain the virtual bonds"
 589       endif
 590       reset_fricmat=1000
 591       if (lang.gt.0) then
 592         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 593         call reada(controlcard,"RWAT",rwat,1.4d0)
 594         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 595         surfarea=index(controlcard,"SURFAREA").gt.0
 596         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 597         if(me.eq.king.or..not.out1file)then
 598          write (iout,'(/a,$)') "Langevin dynamics calculation"
 599          if (lang.eq.1) then
 600           write (iout,'(a/)')
 601      &      " with direct integration of Langevin equations"
 602          else if (lang.eq.2) then
 603           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 604          else if (lang.eq.3) then
 605           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 606          else if (lang.eq.4) then
 607           write (iout,'(a/)') " in overdamped mode"
 608          else
 609           write (iout,'(//a,i5)')
 610      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 611           stop
 612          endif
 613          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 614          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 615          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 616          write (iout,'(a60,f10.5)')
 617      &   "Scaling factor of the friction forces:",scal_fric
 618          if (surfarea) write (iout,'(2a,i10,a)')
 619      &     "Friction coefficients will be scaled by solvent-accessible",
 620      &     " surface area every",reset_fricmat," steps."
 621         endif
 622 c Calculate friction coefficients and bounds of stochastic forces
 623         eta=6*pi*cPoise*etawat
 624         if(me.eq.king.or..not.out1file)
 625      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 626      &   ,eta
 627         gamp=scal_fric*(pstok+rwat)*eta
 628         stdfp=dsqrt(2*Rb*t_bath/d_time)
 629         do i=1,ntyp
 630           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 631           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 632         enddo
 633         if(me.eq.king.or..not.out1file)then
 634          write (iout,'(/2a/)')
 635      &   "Radii of site types and friction coefficients and std's of",
 636      &   " stochastic forces of fully exposed sites"
 637          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 638          do i=1,ntyp
 639           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 640      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 641          enddo
 642         endif
 643       else if (tbf) then
 644         if(me.eq.king.or..not.out1file)then
 645          write (iout,'(a)') "Berendsen bath calculation"
 646          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 647          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 648          if (reset_moment)
 649      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 650      &   count_reset_moment," steps"
 651          if (reset_vel)
 652      &    write (iout,'(a,i10,a)')
 653      &    "Velocities will be reset at random every",count_reset_vel,
 654      &   " steps"
 655         endif
 656       else
 657         if(me.eq.king.or..not.out1file)
 658      &   write (iout,'(a31)') "Microcanonical mode calculation"
 659       endif
 660       if(me.eq.king.or..not.out1file)then
 661        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 662        if (usampl) then
 663           write(iout,*) "MD running with constraints."
 664           write(iout,*) "Equilibration time ", eq_time, " mtus."
 665           write(iout,*) "Constraining ", nfrag," fragments."
 666           write(iout,*) "Length of each fragment, weight and q0:"
 667           do iset=1,nset
 668            write (iout,*) "Set of restraints #",iset
 669            do i=1,nfrag
 670               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 671      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 672            enddo
 673            write(iout,*) "constraints between ", npair, "fragments."
 674            write(iout,*) "constraint pairs, weights and q0:"
 675            do i=1,npair
 676             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 677      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 678            enddo
 679            write(iout,*) "angle constraints within ", nfrag_back,
 680      &      "backbone fragments."
 681            if (loc_qlike) then
 682            write(iout,*) "fragment, weights, q0:"
 683            do i=1,nfrag_back
 684             write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
 685      &         ifrag_back(2,i,iset),
 686      &         wfrag_back(1,i,iset),qin_back(1,i,iset),
 687      &         wfrag_back(2,i,iset),qin_back(2,i,iset),
 688      &         wfrag_back(3,i,iset),qin_back(3,i,iset)
 689            enddo
 690            else
 691            write(iout,*) "fragment, weights:"
 692            do i=1,nfrag_back
 693             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 694      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 695      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 696            enddo
 697            endif
 698           enddo
 699         iset=mod(kolor,nset)+1
 700        endif
 701       endif
 702       if(me.eq.king.or..not.out1file)
 703      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 704       return
 705       end
 706 c------------------------------------------------------------------------------
 707       subroutine molread
 708 C
 709 C Read molecular data.
 710 C
 711       implicit none
 712       include 'DIMENSIONS'
 713 #ifdef MPI
 714       include 'mpif.h'
 715       integer error_msg,ierror,ierr,ierrcode
 716 #endif
 717       include 'COMMON.IOUNITS'
 718       include 'COMMON.GEO'
 719       include 'COMMON.VAR'
 720       include 'COMMON.INTERACT'
 721       include 'COMMON.LOCAL'
 722       include 'COMMON.NAMES'
 723       include 'COMMON.CHAIN'
 724       include 'COMMON.FFIELD'
 725       include 'COMMON.SBRIDGE'
 726       include 'COMMON.HEADER'
 727       include 'COMMON.CONTROL'
 728       include 'COMMON.SAXS'
 729       include 'COMMON.DBASE'
 730       include 'COMMON.THREAD'
 731       include 'COMMON.CONTACTS'
 732       include 'COMMON.TORCNSTR'
 733       include 'COMMON.TIME1'
 734       include 'COMMON.BOUNDS'
 735       include 'COMMON.MD'
 736       include 'COMMON.SETUP'
 737       include 'COMMON.SHIELD'
 738       character*4 sequence(maxres)
 739       integer rescode
 740       double precision x(maxvar)
 741       character*256 pdbfile
 742       character*400 weightcard
 743       character*80 weightcard_t,ucase
 744       integer itype_pdb(maxres)
 745       common /pizda/ itype_pdb
 746       logical seq_comp,fail
 747       double precision energia(0:n_ene)
 748       double precision secprob(3,maxdih_constr)
 749       double precision co
 750       double precision phihel,phibet,sigmahel,sigmabet
 751       integer iti,nsi,maxsi
 752       integer ilen
 753       external ilen
 754       integer iperm,tperm
 755       integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
 756       double precision sumv
 757 C
 758 C Read PDB structure if applicable
 759 C
 760       if (indpdb.gt.0 .or. pdbref) then
 761         read(inp,'(a)') pdbfile
 762         if(me.eq.king.or..not.out1file)
 763      &   write (iout,'(2a)') 'PDB data will be read from file ',
 764      &   pdbfile(:ilen(pdbfile))
 765         open(ipdbin,file=pdbfile,status='old',err=33)
 766         goto 34
 767   33    write (iout,'(a)') 'Error opening PDB file.'
 768         stop
 769   34    continue
 770         call readpdb
 771         do i=1,2*nres
 772           do j=1,3
 773             crefjlee(j,i)=c(j,i)
 774           enddo
 775         enddo
 776 #ifdef DEBUG
 777         do i=1,nres
 778           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 779      &      (crefjlee(j,i+nres),j=1,3)
 780         enddo
 781 #endif
 782         if(me.eq.king.or..not.out1file)
 783      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 784      &   ' nstart_sup=',nstart_sup
 785         do i=1,nres
 786           itype_pdb(i)=itype(i)
 787         enddo
 788         close (ipdbin)
 789         nnt=nstart_sup
 790         nct=nstart_sup+nsup-1
 791         call contact(.false.,ncont_ref,icont_ref,co)
 792
 793         if (sideadd) then
 794          if(me.eq.king.or..not.out1file)
 795      &    write(iout,*)'Adding sidechains'
 796          maxsi=1000
 797          do i=2,nres-1
 798           iti=itype(i)
 799           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 800             nsi=0
 801             fail=.true.
 802             do while (fail.and.nsi.le.maxsi)
 803               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 804               nsi=nsi+1
 805             enddo
 806             if(fail) write(iout,*)'Adding sidechain failed for res ',
 807      &              i,' after ',nsi,' trials'
 808           endif
 809          enddo
 810         endif
 811       endif
 812       if (indpdb.eq.0) then
 813 C Read sequence if not taken from the pdb file.
 814         read (inp,*) nres
 815 c        print *,'nres=',nres
 816         if (iscode.gt.0) then
 817           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 818         else
 819           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 820         endif
 821 C Convert sequence to numeric code
 822         do i=1,nres
 823           itype(i)=rescode(i,sequence(i),iscode)
 824         enddo
 825       endif
 826 c      print *,nres
 827 c      print '(20i4)',(itype(i),i=1,nres)
 828       do i=1,nres
 829 #ifdef PROCOR
 830         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 831 #else
 832         if (itype(i).eq.ntyp1) then
 833 #endif
 834           itel(i)=0
 835 #ifdef PROCOR
 836         else if (iabs(itype(i+1)).ne.20) then
 837 #else
 838         else if (iabs(itype(i)).ne.20) then
 839 #endif
 840           itel(i)=1
 841         else
 842           itel(i)=2
 843         endif
 844       enddo
 845       if(me.eq.king.or..not.out1file)then
 846        write (iout,*) "ITEL"
 847        do i=1,nres-1
 848          write (iout,*) i,itype(i),itel(i)
 849        enddo
 850 c       print *,'Call Read_Bridge.'
 851       endif
 852       nnt=1
 853       nct=nres
 854 cd      print *,'NNT=',NNT,' NCT=',NCT
 855       call seq2chains(nres,itype,nchain,chain_length,chain_border,
 856      &  ireschain)
 857       chain_border1(1,1)=1
 858       chain_border1(2,1)=chain_border(2,1)+1
 859       do i=2,nchain-1
 860         chain_border1(1,i)=chain_border(1,i)-1
 861         chain_border1(2,i)=chain_border(2,i)+1
 862       enddo
 863       if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
 864       chain_border1(2,nchain)=nres
 865       write(iout,*) "nres",nres," nchain",nchain
 866       do i=1,nchain
 867         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
 868      &    chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
 869       enddo
 870       call chain_symmetry(nchain,nres,itype,chain_border,
 871      &    chain_length,npermchain,tabpermchain,nchain_group,nequiv,
 872      &    iequiv,mapchain)
 873 c      do i=1,nres
 874 c        write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 875 c     &    ii=1,npermchain)
 876 c      enddo
 877 #ifdef DEBUG
 878       write(iout,*) "residue permutations"
 879       do i=1,nres
 880         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
 881       enddo
 882       call flush(iout)
 883 #endif
 884       if (itype(1).eq.ntyp1) nnt=2
 885       if (itype(nres).eq.ntyp1) nct=nct-1
 886       write (iout,*) "nnt",nnt," nct",nct
 887       call flush(iout)
 888 #ifdef DFA
 889       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 890      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 891        call init_dfa_vars
 892 c       print*, 'init_dfa_vars finished!'
 893        call read_dfa_info
 894 c       print*, 'read_dfa_info finished!'
 895       endif
 896 #endif
 897       if (pdbref) then
 898         if(me.eq.king.or..not.out1file)
 899      &   write (iout,'(a,i3)') 'nsup=',nsup
 900         nstart_seq=nnt
 901         if (nsup.le.(nct-nnt+1)) then
 902           do i=0,nct-nnt+1-nsup
 903             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 904               nstart_seq=nnt+i
 905               goto 111
 906             endif
 907           enddo
 908           write (iout,'(a)')
 909      &            'Error - sequences to be superposed do not match.'
 910           stop
 911         else
 912           do i=0,nsup-(nct-nnt+1)
 913             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 914      &      then
 915               nstart_sup=nstart_sup+i
 916               nsup=nct-nnt+1
 917               goto 111
 918             endif
 919           enddo
 920           write (iout,'(a)')
 921      &            'Error - sequences to be superposed do not match.'
 922         endif
 923   111   continue
 924         if (nsup.eq.0) nsup=nct-nnt
 925         if (nstart_sup.eq.0) nstart_sup=nnt
 926         if (nstart_seq.eq.0) nstart_seq=nnt
 927         if(me.eq.king.or..not.out1file)
 928      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 929      &                 ' nstart_seq=',nstart_seq
 930       endif
 931       call read_bridge
 932 C 8/13/98 Set limits to generating the dihedral angles
 933       do i=1,nres
 934         phibound(1,i)=-pi
 935         phibound(2,i)=pi
 936       enddo
 937       read (inp,*) ndih_constr
 938       write (iout,*) "ndih_constr",ndih_constr
 939       if (ndih_constr.gt.0) then
 940         raw_psipred=.false.
 941 C        read (inp,*) ftors
 942         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 943      &  i=1,ndih_constr)
 944 #ifdef MPI
 945         if(me.eq.king.or..not.out1file)then
 946 #endif
 947          write (iout,*)
 948      &   'There are',ndih_constr,' restraints on gamma angles.'
 949          do i=1,ndih_constr
 950           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 951      &    ftors(i)
 952          enddo
 953 #ifdef MPI
 954         endif
 955 #endif
 956         do i=1,ndih_constr
 957           phi0(i)=deg2rad*phi0(i)
 958           drange(i)=deg2rad*drange(i)
 959         enddo
 960 C        if(me.eq.king.or..not.out1file)
 961 C     &   write (iout,*) 'FTORS',ftors
 962         do i=1,ndih_constr
 963           ii = idih_constr(i)
 964           phibound(1,ii) = phi0(i)-drange(i)
 965           phibound(2,ii) = phi0(i)+drange(i)
 966         enddo
 967 #ifdef MPI
 968         if (me.eq.king .or. .not.out1file) then
 969 #endif
 970          write (iout,'(a)') 'Boundaries in gamma angle sampling:'
 971          do i=1,nres
 972            write (iout,'(a3,i5,2f10.1)')
 973      &    restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 974          enddo
 975 #ifdef MP
 976         endif
 977 #endif
 978       else if (ndih_constr.lt.0) then
 979         raw_psipred=.true.
 980         call card_concat(weightcard)
 981         call reada(weightcard,"PHIHEL",phihel,50.0D0)
 982         call reada(weightcard,"PHIBET",phibet,180.0D0)
 983         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
 984         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
 985         call reada(weightcard,"WDIHC",wdihc,0.591D0)
 986         write (iout,*) "Weight of dihedral angle restraints",wdihc
 987         read(inp,'(9x,3f7.3)')
 988 c     &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 989      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 990         write (iout,*) "The secprob array"
 991         do i=nnt,nct
 992           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
 993         enddo
 994         ndih_constr=0
 995         iconstr_dih=0
 996         do i=nnt+3,nct
 997           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
 998      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
 999             ndih_constr=ndih_constr+1
1000             idih_constr(ndih_constr)=i
1001             iconstr_dih(i)=ndih_constr
1002             sumv=0.0d0
1003             do j=1,3
1004               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
1005               sumv=sumv+vpsipred(j,ndih_constr)
1006             enddo
1007             if (sumv.gt.0.0d0) then
1008               do j=1,3
1009                 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1010               enddo
1011             else
1012               vpsipred(1,ndih_constr)=1.0d0
1013               vpsipred(2,ndih_constr)=0.0d0
1014               vpsipred(3,ndih_constr)=0.0d0
1015             endif
1016             phibound(1,ndih_constr)=phihel*deg2rad
1017             phibound(2,ndih_constr)=phibet*deg2rad
1018             sdihed(1,ndih_constr)=sigmahel*deg2rad
1019             sdihed(2,ndih_constr)=sigmabet*deg2rad
1020           endif
1021         enddo
1022 #ifdef MPI
1023         if(me.eq.king.or..not.out1file)then
1024 #endif
1025          write (iout,*)
1026      &   'There are',ndih_constr,
1027      &   ' bimodal restraints on gamma angles.'
1028          do i=1,ndih_constr
1029            write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1030      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1031      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1032      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1033      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1034      &      (vpsipred(j,i),j=1,3)
1035          enddo
1036 #ifdef MPI
1037         endif
1038 #endif
1039 c Raw psipred input
1040       endif
1041 C first setting the theta boundaries to 0 to pi
1042 C this mean that there is no energy penalty for any angle occuring this can be applied
1043 C for generate random conformation but is not implemented in this way
1044 C      do i=1,nres
1045 C        thetabound(1,i)=0
1046 C        thetabound(2,i)=pi
1047 C      enddo
1048 C begin reading theta constrains this is quartic constrains allowing to
1049 C have smooth second derivative
1050       if (with_theta_constr) then
1051 C with_theta_constr is keyword allowing for occurance of theta constrains
1052       read (inp,*) ntheta_constr
1053 C ntheta_constr is the number of theta constrains
1054       if (ntheta_constr.gt.0) then
1055 C        read (inp,*) ftors
1056         read (inp,*) (itheta_constr(i),theta_constr0(i),
1057      &  theta_drange(i),for_thet_constr(i),
1058      &  i=1,ntheta_constr)
1059 C the above code reads from 1 to ntheta_constr
1060 C itheta_constr(i) residue i for which is theta_constr
1061 C theta_constr0 the global minimum value
1062 C theta_drange is range for which there is no energy penalty
1063 C for_thet_constr is the force constant for quartic energy penalty
1064 C E=k*x**4
1065         if(me.eq.king.or..not.out1file)then
1066          write (iout,*)
1067      &   'There are',ntheta_constr,' constraints on phi angles.'
1068          do i=1,ntheta_constr
1069           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1070      &    theta_drange(i),
1071      &    for_thet_constr(i)
1072          enddo
1073         endif
1074         do i=1,ntheta_constr
1075           theta_constr0(i)=deg2rad*theta_constr0(i)
1076           theta_drange(i)=deg2rad*theta_drange(i)
1077         enddo
1078 C        if(me.eq.king.or..not.out1file)
1079 C     &   write (iout,*) 'FTORS',ftors
1080 C        do i=1,ntheta_constr
1081 C          ii = itheta_constr(i)
1082 C          thetabound(1,ii) = phi0(i)-drange(i)
1083 C          thetabound(2,ii) = phi0(i)+drange(i)
1084 C        enddo
1085       endif ! ntheta_constr.gt.0
1086       endif! with_theta_constr
1087 C
1088 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1089 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
1090 c--- Zscore rms -------
1091       if (nz_start.eq.0) nz_start=nnt
1092       if (nz_end.eq.0 .and. nsup.gt.0) then
1093         nz_end=nnt+nsup-1
1094       else if (nz_end.eq.0) then
1095         nz_end=nct
1096       endif
1097       if(me.eq.king.or..not.out1file)then
1098        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1099        write (iout,*) 'IZ_SC=',iz_sc
1100       endif
1101 c----------------------
1102       call init_int_table
1103       if (refstr) then
1104         if (.not.pdbref) then
1105           call read_angles(inp,*38)
1106           call bond_regular
1107           goto 39
1108    38     write (iout,'(a)') 'Error reading reference structure.'
1109 #ifdef MPI
1110           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1111           stop 'Error reading reference structure'
1112 #endif
1113    39     call chainbuild_extconf
1114           call setup_var
1115 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1116           nstart_sup=nnt
1117           nstart_seq=nnt
1118           nsup=nct-nnt+1
1119           do i=1,2*nres
1120             do j=1,3
1121               cref(j,i)=c(j,i)
1122             enddo
1123           enddo
1124           call contact(.true.,ncont_ref,icont_ref,co)
1125         endif
1126       endif
1127 c        print *, "A TU"
1128 c      write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1129       call flush(iout)
1130       if (constr_dist.gt.0) call read_dist_constr
1131 c      write (iout,*) "After read_dist_constr nhpb",nhpb
1132       if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1133       call hpb_partition
1134       call NMRpeak_partition
1135       if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1136       if (pdbref) then
1137         if(me.eq.king.or..not.out1file)
1138      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1139         do i=1,ncont_ref
1140           do j=1,2
1141             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1142           enddo
1143           if(me.eq.king.or..not.out1file)
1144      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1145      &     icont_ref(1,i),' ',
1146      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1147         enddo
1148       endif
1149       write (iout,*) "calling read_saxs_consrtr",nsaxs
1150       if (nsaxs.gt.0) call read_saxs_constr
1151 c      write (iout,*) "After read_saxs_constr"
1152 c      call flush(iout)
1153       if (constr_homology.gt.0) then
1154 c        write (iout,*) "Calling read_constr_homology"
1155 c        call flush(iout)
1156         call read_constr_homology
1157         if (indpdb.gt.0 .or. pdbref) then
1158           do i=1,2*nres
1159             do j=1,3
1160               c(j,i)=crefjlee(j,i)
1161               cref(j,i)=crefjlee(j,i)
1162             enddo
1163           enddo
1164         endif
1165 #ifdef DEBUG
1166         write (iout,*) "sc_loc_geom: Array C"
1167         do i=1,nres
1168           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1169      &      (c(j,i+nres),j=1,3)
1170         enddo
1171         write (iout,*) "Array Cref"
1172         do i=1,nres
1173           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1174      &      (cref(j,i+nres),j=1,3)
1175         enddo
1176 #endif
1177        call int_from_cart1(.false.)
1178        call sc_loc_geom(.false.)
1179        do i=1,nres
1180          thetaref(i)=theta(i)
1181          phiref(i)=phi(i)
1182        enddo
1183        do i=1,nres-1
1184          do j=1,3
1185            dc(j,i)=c(j,i+1)-c(j,i)
1186            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1187          enddo
1188        enddo
1189        do i=2,nres-1
1190          do j=1,3
1191            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1192            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1193          enddo
1194        enddo
1195       else
1196         homol_nset=0
1197         if (start_from_model) then
1198           nmodel_start=0
1199           do
1200             read(inp,'(a)',end=332,err=332) pdbfile
1201             if (me.eq.king .or. .not. out1file)
1202      &        write (iout,'(a,5x,a)') 'Opening PDB file',
1203      &        pdbfile(:ilen(pdbfile))
1204             open(ipdbin,file=pdbfile,status='old',err=336)
1205             goto 335
1206  336        write (iout,'(a,5x,a)') 'Error opening PDB file',
1207      &      pdbfile(:ilen(pdbfile))
1208             call flush(iout)
1209             stop
1210  335        continue
1211             unres_pdb=.false.
1212             nres_temp=nres
1213 c            call readpdb
1214             call readpdb_template(nmodel_start+1)
1215             close(ipdbin)
1216             if (nres.ge.nres_temp) then
1217               nmodel_start=nmodel_start+1
1218               pdbfiles_chomo(nmodel_start)=pdbfile
1219               do i=1,2*nres
1220                 do j=1,3
1221                   chomo(j,i,nmodel_start)=c(j,i)
1222                 enddo
1223               enddo
1224             else
1225               if (me.eq.king .or. .not. out1file)
1226      &          write (iout,'(a,2i7,1x,a)')
1227      &           "Different number of residues",nres_temp,nres,
1228      &           " model skipped."
1229             endif
1230             nres=nres_temp
1231           enddo
1232   332     continue
1233           if (nmodel_start.eq.0) then
1234             if (me.eq.king .or. .not. out1file)
1235      &        write (iout,'(a)')
1236      &        "No valid starting model found START_FROM_MODELS is OFF"
1237               start_from_model=.false.
1238           endif
1239           write (iout,*) "nmodel_start",nmodel_start
1240         endif
1241       endif
1242
1243
1244 C      endif
1245       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1246      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1247      &    modecalc.ne.10) then
1248 C If input structure hasn't been supplied from the PDB file read or generate
1249 C initial geometry.
1250         if (iranconf.eq.0 .and. .not. extconf .and. .not.
1251      &     start_from_model) then
1252           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1253      &     write (iout,'(a)') 'Initial geometry will be read in.'
1254           if (read_cart) then
1255             read(inp,'(8f10.5)',end=36,err=36)
1256      &       ((c(l,k),l=1,3),k=1,nres),
1257      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1258             if (nnt.gt.1) c(:,nres+1)=c(:,1)
1259             if (nct.lt.nres) c(:,2*nres)=c(:,nres)
1260 c            write (iout,*) "Exit READ_CART"
1261 c            write (iout,'(8f10.5)')
1262 c     &       ((c(l,k),l=1,3),k=1,nres),
1263 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1264             call cartprint
1265             do j=1,3
1266               dc(j,0)=c(j,1)
1267             enddo
1268             do i=1,nres-1
1269               do j=1,3
1270                 dc(j,i)=c(j,i+1)-c(j,i)
1271 c                dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1272               enddo
1273             enddo
1274             do i=nnt,nct
1275               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1276                 do j=1,3
1277                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1278 c                  dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1279                 enddo
1280               else
1281                 do j=1,3
1282                   dc(j,i+nres)=0.0d0
1283 c                  dc_norm(j,i+nres)=0.0d0
1284                 enddo
1285               endif
1286             enddo
1287             call int_from_cart1(.true.)
1288             write (iout,*) "Finish INT_TO_CART"
1289 c      write (iout,*) "DC"
1290 c      do i=1,nres
1291 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1292 c     &    (dc(j,i+nres),j=1,3)
1293 c      enddo
1294 c            call cartprint
1295 c            call setup_var
1296 c            return
1297           else
1298             call read_angles(inp,*36)
1299             call bond_regular
1300             call chainbuild_extconf
1301           endif
1302           goto 37
1303    36     write (iout,'(a)') 'Error reading angle file.'
1304 #ifdef MPI
1305           call mpi_finalize( MPI_COMM_WORLD,IERR )
1306 #endif
1307           stop 'Error reading angle file.'
1308    37     continue
1309         else if (extconf) then
1310           if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1311      &      write (iout,'(a)') 'Extended chain initial geometry.'
1312           do i=3,nres
1313             theta(i)=90d0*deg2rad
1314           enddo
1315           do i=4,nres
1316             phi(i)=180d0*deg2rad
1317           enddo
1318           do i=2,nres-1
1319             alph(i)=110d0*deg2rad
1320           enddo
1321           do i=2,nres-1
1322             omeg(i)=-120d0*deg2rad
1323             if (itype(i).le.0) omeg(i)=-omeg(i)
1324           enddo
1325           call bond_regular
1326           call chainbuild_extconf
1327         else if (.not. start_from_model) then
1328           if(me.eq.king.or..not.out1file)
1329      &     write (iout,'(a)') 'Random-generated initial geometry.'
1330           call bond_regular
1331 #ifdef MPI
1332           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1333      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1334 #endif
1335             do itrial=1,100
1336               itmp=1
1337               call gen_rand_conf(itmp,*30)
1338               goto 40
1339    30         write (iout,*) 'Failed to generate random conformation',
1340      &          ', itrial=',itrial
1341               write (*,*) 'Processor:',me,
1342      &          ' Failed to generate random conformation',
1343      &          ' itrial=',itrial
1344               call intout
1345
1346 #ifdef AIX
1347               call flush_(iout)
1348 #else
1349               call flush(iout)
1350 #endif
1351             enddo
1352             write (iout,'(a,i3,a)') 'Processor:',me,
1353      &        ' error in generating random conformation.'
1354             write (*,'(a,i3,a)') 'Processor:',me,
1355      &        ' error in generating random conformation.'
1356             call flush(iout)
1357 #ifdef MPI
1358             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1359    40       continue
1360           endif
1361 #else
1362           write (*,'(a)')
1363      &      ' error in generating random conformation.'
1364           stop
1365    40     continue
1366 #endif
1367         endif
1368       elseif (modecalc.eq.4) then
1369         read (inp,'(a)') intinname
1370         open (intin,file=intinname,status='old',err=333)
1371         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1372      &  write (iout,'(a)') 'intinname',intinname
1373         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1374         goto 334
1375   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1376 #ifdef MPI
1377         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1378 #endif
1379         stop 'Error opening angle file.'
1380   334   continue
1381
1382       endif
1383 C Generate distance constraints, if the PDB structure is to be regularized.
1384       if (nthread.gt.0) then
1385         call read_threadbase
1386       endif
1387       call setup_var
1388       if (me.eq.king .or. .not. out1file)
1389      & call intout
1390       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1391         write (iout,'(/a,i3,a)')
1392      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1393         write (iout,'(20i4)') (iss(i),i=1,ns)
1394        if (dyn_ss) then
1395           write(iout,*)"Running with dynamic disulfide-bond formation"
1396        else
1397         write (iout,'(/a/)') 'Pre-formed links are:'
1398         do i=1,nss
1399           i1=ihpb(i)-nres
1400           i2=jhpb(i)-nres
1401           it1=itype(i1)
1402           it2=itype(i2)
1403           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1404      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1405      &    ebr,forcon(i)
1406         enddo
1407         write (iout,'(a)')
1408        endif
1409       endif
1410       if (ns.gt.0.and.dyn_ss) then
1411           do i=nss+1,nhpb
1412             ihpb(i-nss)=ihpb(i)
1413             jhpb(i-nss)=jhpb(i)
1414             forcon(i-nss)=forcon(i)
1415             dhpb(i-nss)=dhpb(i)
1416           enddo
1417           nhpb=nhpb-nss
1418           nss=0
1419           call hpb_partition
1420           do i=1,ns
1421             dyn_ss_mask(iss(i))=.true.
1422           enddo
1423       endif
1424 c --- test
1425 c      call cartprint
1426 c      write (iout,*) "DC"
1427 c      do i=1,nres
1428 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1429 c     &    (dc(j,i+nres),j=1,3)
1430 c      enddo
1431       if (i2ndstr.gt.0) call secstrp2dihc
1432 c      call geom_to_var(nvar,x)
1433 c      call etotal(energia(0))
1434 c      call enerprint(energia(0))
1435 c      call briefout(0,etot)
1436 c      stop
1437 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1438 cd    write (iout,'(a)') 'Variable list:'
1439 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1440 #ifdef MPI
1441       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1442      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1443      &  'Processor',myrank,': end reading molecular data.'
1444 #endif
1445 c      print *,"A TU?"
1446       return
1447       end
1448 c--------------------------------------------------------------------------
1449       logical function seq_comp(itypea,itypeb,length)
1450       implicit none
1451       integer length,itypea(length),itypeb(length)
1452       integer i
1453       do i=1,length
1454         if (itypea(i).ne.itypeb(i)) then
1455           seq_comp=.false.
1456           return
1457         endif
1458       enddo
1459       seq_comp=.true.
1460       return
1461       end
1462 c-----------------------------------------------------------------------------
1463       subroutine read_bridge
1464 C Read information about disulfide bridges.
1465       implicit none
1466       include 'DIMENSIONS'
1467 #ifdef MPI
1468       include 'mpif.h'
1469       integer ierror
1470 #endif
1471       include 'COMMON.IOUNITS'
1472       include 'COMMON.GEO'
1473       include 'COMMON.VAR'
1474       include 'COMMON.INTERACT'
1475       include 'COMMON.LOCAL'
1476       include 'COMMON.NAMES'
1477       include 'COMMON.CHAIN'
1478       include 'COMMON.FFIELD'
1479       include 'COMMON.SBRIDGE'
1480       include 'COMMON.HEADER'
1481       include 'COMMON.CONTROL'
1482       include 'COMMON.DBASE'
1483       include 'COMMON.THREAD'
1484       include 'COMMON.TIME1'
1485       include 'COMMON.SETUP'
1486       integer i,j
1487 C Read bridging residues.
1488       read (inp,*) ns,(iss(i),i=1,ns)
1489 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1490 c     numbers
1491       icys=0
1492       do i=1,ns
1493         icys(iss(i))=i
1494       enddo
1495 c      print *,'ns=',ns
1496       if(me.eq.king.or..not.out1file)
1497      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1498 C Check whether the specified bridging residues are cystines.
1499       do i=1,ns
1500         if (itype(iss(i)).ne.1) then
1501           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1502      &   'Do you REALLY think that the residue ',
1503      &    restyp(itype(iss(i))),i,
1504      &   ' can form a disulfide bridge?!!!'
1505           write (*,'(2a,i3,a)')
1506      &   'Do you REALLY think that the residue ',
1507      &    restyp(itype(iss(i))),i,
1508      &   ' can form a disulfide bridge?!!!'
1509 #ifdef MPI
1510          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1511          stop
1512 #endif
1513         endif
1514       enddo
1515 C Read preformed bridges.
1516       if (ns.gt.0) then
1517       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1518       if(fg_rank.eq.0)
1519      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1520       if (nss.gt.0) then
1521         nhpb=nss
1522 C Check if the residues involved in bridges are in the specified list of
1523 C bridging residues.
1524         do i=1,nss
1525           do j=1,i-1
1526             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1527      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1528               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1529      &      ' contains residues present in other pairs.'
1530               write (*,'(a,i3,a)') 'Disulfide pair',i,
1531      &      ' contains residues present in other pairs.'
1532 #ifdef MPI
1533               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1534               stop
1535 #endif
1536             endif
1537           enddo
1538           do j=1,ns
1539             if (ihpb(i).eq.iss(j)) goto 10
1540           enddo
1541           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1542    10     continue
1543           do j=1,ns
1544             if (jhpb(i).eq.iss(j)) goto 20
1545           enddo
1546           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1547    20     continue
1548 c          dhpb(i)=dbr
1549 c          forcon(i)=fbr
1550         enddo
1551         do i=1,nss
1552           ihpb(i)=ihpb(i)+nres
1553           jhpb(i)=jhpb(i)+nres
1554         enddo
1555       endif
1556       endif
1557       return
1558       end
1559 c----------------------------------------------------------------------------
1560       subroutine read_x(kanal,*)
1561       implicit none
1562       include 'DIMENSIONS'
1563       include 'COMMON.GEO'
1564       include 'COMMON.VAR'
1565       include 'COMMON.CHAIN'
1566       include 'COMMON.IOUNITS'
1567       include 'COMMON.CONTROL'
1568       include 'COMMON.LOCAL'
1569       include 'COMMON.INTERACT'
1570       integer i,j,k,l,kanal
1571 c Read coordinates from input
1572 c
1573       read(kanal,'(8f10.5)',end=10,err=10)
1574      &  ((c(l,k),l=1,3),k=1,nres),
1575      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1576       do j=1,3
1577         c(j,nres+1)=c(j,1)
1578         c(j,2*nres)=c(j,nres)
1579       enddo
1580       call int_from_cart1(.false.)
1581       do i=1,nres-1
1582         do j=1,3
1583           dc(j,i)=c(j,i+1)-c(j,i)
1584           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1585         enddo
1586       enddo
1587       do i=nnt,nct
1588         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1589           do j=1,3
1590             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1591             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1592           enddo
1593         endif
1594       enddo
1595
1596       return
1597    10 return1
1598       end
1599 c----------------------------------------------------------------------------
1600       subroutine read_threadbase
1601       implicit none
1602       include 'DIMENSIONS'
1603       include 'COMMON.IOUNITS'
1604       include 'COMMON.GEO'
1605       include 'COMMON.VAR'
1606       include 'COMMON.INTERACT'
1607       include 'COMMON.LOCAL'
1608       include 'COMMON.NAMES'
1609       include 'COMMON.CHAIN'
1610       include 'COMMON.FFIELD'
1611       include 'COMMON.SBRIDGE'
1612       include 'COMMON.HEADER'
1613       include 'COMMON.CONTROL'
1614       include 'COMMON.DBASE'
1615       include 'COMMON.THREAD'
1616       include 'COMMON.TIME1'
1617       integer i,j,k
1618       double precision dist
1619 C Read pattern database for threading.
1620       read (icbase,*) nseq
1621       do i=1,nseq
1622         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1623      &   nres_base(2,i),nres_base(3,i)
1624         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1625      &   nres_base(1,i))
1626 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1627 c    &   nres_base(2,i),nres_base(3,i)
1628 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1629 c    &   nres_base(1,i))
1630       enddo
1631       close (icbase)
1632       if (weidis.eq.0.0D0) weidis=0.1D0
1633       do i=nnt,nct
1634         do j=i+2,nct
1635           nhpb=nhpb+1
1636           ihpb(nhpb)=i
1637           jhpb(nhpb)=j
1638           forcon(nhpb)=weidis
1639         enddo
1640       enddo
1641       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1642       write (iout,'(a,i5)') 'nexcl: ',nexcl
1643       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1644       return
1645       end
1646 c------------------------------------------------------------------------------
1647       subroutine setup_var
1648       implicit none
1649       include 'DIMENSIONS'
1650       include 'COMMON.IOUNITS'
1651       include 'COMMON.GEO'
1652       include 'COMMON.VAR'
1653       include 'COMMON.INTERACT'
1654       include 'COMMON.LOCAL'
1655       include 'COMMON.NAMES'
1656       include 'COMMON.CHAIN'
1657       include 'COMMON.FFIELD'
1658       include 'COMMON.SBRIDGE'
1659       include 'COMMON.HEADER'
1660       include 'COMMON.CONTROL'
1661       include 'COMMON.DBASE'
1662       include 'COMMON.THREAD'
1663       include 'COMMON.TIME1'
1664       integer i
1665 C Set up variable list.
1666       ntheta=nres-2
1667       nphi=nres-3
1668       nvar=ntheta+nphi
1669       nside=0
1670       do i=2,nres-1
1671         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1672           nside=nside+1
1673           ialph(i,1)=nvar+nside
1674           ialph(nside,2)=i
1675         endif
1676       enddo
1677       if (indphi.gt.0) then
1678         nvar=nphi
1679       else if (indback.gt.0) then
1680         nvar=nphi+ntheta
1681       else
1682         nvar=nvar+2*nside
1683       endif
1684       return
1685       end
1686 c----------------------------------------------------------------------------
1687       subroutine gen_dist_constr
1688 C Generate CA distance constraints.
1689       implicit none
1690       include 'DIMENSIONS'
1691       include 'COMMON.IOUNITS'
1692       include 'COMMON.GEO'
1693       include 'COMMON.VAR'
1694       include 'COMMON.INTERACT'
1695       include 'COMMON.LOCAL'
1696       include 'COMMON.NAMES'
1697       include 'COMMON.CHAIN'
1698       include 'COMMON.FFIELD'
1699       include 'COMMON.SBRIDGE'
1700       include 'COMMON.HEADER'
1701       include 'COMMON.CONTROL'
1702       include 'COMMON.DBASE'
1703       include 'COMMON.THREAD'
1704       include 'COMMON.SPLITELE'
1705       include 'COMMON.TIME1'
1706       integer i,j,itype_pdb(maxres)
1707       common /pizda/ itype_pdb
1708       double precision dd
1709       double precision dist
1710       character*2 iden
1711 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1712 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1713 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1714 cd     & ' nsup',nsup
1715       do i=nstart_sup,nstart_sup+nsup-1
1716 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1717 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1718         do j=i+2,nstart_sup+nsup-1
1719 c 5/24/2020 Adam: Cutoff included to reduce array size
1720           dd = dist(i,j)
1721           if (dd.gt.r_cut_int) cycle
1722           nhpb=nhpb+1
1723           ihpb(nhpb)=i+nstart_seq-nstart_sup
1724           jhpb(nhpb)=j+nstart_seq-nstart_sup
1725           forcon(nhpb)=weidis
1726           dhpb(nhpb)=dd
1727         enddo
1728       enddo
1729 cd      write (iout,'(a)') 'Distance constraints:'
1730 cd      do i=nss+1,nhpb
1731 cd        ii=ihpb(i)
1732 cd        jj=jhpb(i)
1733 cd        iden='CA'
1734 cd        if (ii.gt.nres) then
1735 cd          iden='SC'
1736 cd          ii=ii-nres
1737 cd          jj=jj-nres
1738 cd        endif
1739 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1740 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1741 cd     &  dhpb(i),forcon(i)
1742 cd      enddo
1743       return
1744       end
1745 c----------------------------------------------------------------------------
1746       subroutine map_read
1747       implicit none
1748       include 'DIMENSIONS'
1749       include 'COMMON.MAP'
1750       include 'COMMON.IOUNITS'
1751       integer imap
1752       character*3 angid(4) /'THE','PHI','ALP','OME'/
1753       character*80 mapcard,ucase
1754       do imap=1,nmap
1755         read (inp,'(a)') mapcard
1756         mapcard=ucase(mapcard)
1757         if (index(mapcard,'PHI').gt.0) then
1758           kang(imap)=1
1759         else if (index(mapcard,'THE').gt.0) then
1760           kang(imap)=2
1761         else if (index(mapcard,'ALP').gt.0) then
1762           kang(imap)=3
1763         else if (index(mapcard,'OME').gt.0) then
1764           kang(imap)=4
1765         else
1766           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1767           stop 'Error - illegal variable spec in MAP card.'
1768         endif
1769         call readi (mapcard,'RES1',res1(imap),0)
1770         call readi (mapcard,'RES2',res2(imap),0)
1771         if (res1(imap).eq.0) then
1772           res1(imap)=res2(imap)
1773         else if (res2(imap).eq.0) then
1774           res2(imap)=res1(imap)
1775         endif
1776         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1777           write (iout,'(a)')
1778      &    'Error - illegal definition of variable group in MAP.'
1779           stop 'Error - illegal definition of variable group in MAP.'
1780         endif
1781         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1782         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1783         call readi(mapcard,'NSTEP',nstep(imap),0)
1784         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1785           write (iout,'(a)')
1786      &     'Illegal boundary and/or step size specification in MAP.'
1787           stop 'Illegal boundary and/or step size specification in MAP.'
1788         endif
1789       enddo ! imap
1790       return
1791       end
1792 c----------------------------------------------------------------------------
1793       subroutine csaread
1794       implicit none
1795       include 'DIMENSIONS'
1796       include 'COMMON.IOUNITS'
1797       include 'COMMON.GEO'
1798       include 'COMMON.CSA'
1799       include 'COMMON.BANK'
1800       include 'COMMON.CONTROL'
1801       character*80 ucase
1802       character*620 mcmcard
1803       call card_concat(mcmcard)
1804
1805       call readi(mcmcard,'NCONF',nconf,50)
1806       call readi(mcmcard,'NADD',nadd,0)
1807       call readi(mcmcard,'JSTART',jstart,1)
1808       call readi(mcmcard,'JEND',jend,1)
1809       call readi(mcmcard,'NSTMAX',nstmax,500000)
1810       call readi(mcmcard,'N0',n0,1)
1811       call readi(mcmcard,'N1',n1,6)
1812       call readi(mcmcard,'N2',n2,4)
1813       call readi(mcmcard,'N3',n3,0)
1814       call readi(mcmcard,'N4',n4,0)
1815       call readi(mcmcard,'N5',n5,0)
1816       call readi(mcmcard,'N6',n6,10)
1817       call readi(mcmcard,'N7',n7,0)
1818       call readi(mcmcard,'N8',n8,0)
1819       call readi(mcmcard,'N9',n9,0)
1820       call readi(mcmcard,'N14',n14,0)
1821       call readi(mcmcard,'N15',n15,0)
1822       call readi(mcmcard,'N16',n16,0)
1823       call readi(mcmcard,'N17',n17,0)
1824       call readi(mcmcard,'N18',n18,0)
1825
1826       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1827
1828       call readi(mcmcard,'NDIFF',ndiff,2)
1829       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1830       call readi(mcmcard,'IS1',is1,1)
1831       call readi(mcmcard,'IS2',is2,8)
1832       call readi(mcmcard,'NRAN0',nran0,4)
1833       call readi(mcmcard,'NRAN1',nran1,2)
1834       call readi(mcmcard,'IRR',irr,1)
1835       call readi(mcmcard,'NSEED',nseed,20)
1836       call readi(mcmcard,'NTOTAL',ntotal,10000)
1837       call reada(mcmcard,'CUT1',cut1,2.0d0)
1838       call reada(mcmcard,'CUT2',cut2,5.0d0)
1839       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1840       call readi(mcmcard,'ICMAX',icmax,3)
1841       call readi(mcmcard,'IRESTART',irestart,0)
1842 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1843       ntbankm=0
1844 c!bankt
1845       call reada(mcmcard,'DELE',dele,20.0d0)
1846       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1847       call readi(mcmcard,'IREF',iref,0)
1848       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1849       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1850       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1851       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1852       write (iout,*) "NCONF_IN",nconf_in
1853       return
1854       end
1855 c----------------------------------------------------------------------------
1856       subroutine mcmread
1857       implicit none
1858       include 'DIMENSIONS'
1859       include 'COMMON.MCM'
1860       include 'COMMON.MCE'
1861       include 'COMMON.IOUNITS'
1862       character*80 ucase
1863       character*320 mcmcard
1864       integer i
1865       call card_concat(mcmcard)
1866       call readi(mcmcard,'MAXACC',maxacc,100)
1867       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1868       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1869       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1870       call readi(mcmcard,'MAXREPM',maxrepm,200)
1871       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1872       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1873       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1874       call reada(mcmcard,'E_UP',e_up,5.0D0)
1875       call reada(mcmcard,'DELTE',delte,0.1D0)
1876       call readi(mcmcard,'NSWEEP',nsweep,5)
1877       call readi(mcmcard,'NSTEPH',nsteph,0)
1878       call readi(mcmcard,'NSTEPC',nstepc,0)
1879       call reada(mcmcard,'TMIN',tmin,298.0D0)
1880       call reada(mcmcard,'TMAX',tmax,298.0D0)
1881       call readi(mcmcard,'NWINDOW',nwindow,0)
1882       call readi(mcmcard,'PRINT_MC',print_mc,0)
1883       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1884       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1885       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1886       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1887       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1888       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1889       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1890       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1891       if (nwindow.gt.0) then
1892         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1893         do i=1,nwindow
1894           winlen(i)=winend(i)-winstart(i)+1
1895         enddo
1896       endif
1897       if (tmax.lt.tmin) tmax=tmin
1898       if (tmax.eq.tmin) then
1899         nstepc=0
1900         nsteph=0
1901       endif
1902       if (nstepc.gt.0 .and. nsteph.gt.0) then
1903         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1904         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1905       endif
1906 C Probabilities of different move types
1907       sumpro_type(0)=0.0D0
1908       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1909       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1910       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1911       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1912       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1913       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1914       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1915       do i=1,MaxMoveType
1916         print *,'i',i,' sumprotype',sumpro_type(i)
1917         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1918         print *,'i',i,' sumprotype',sumpro_type(i)
1919       enddo
1920       return
1921       end
1922 c----------------------------------------------------------------------------
1923       subroutine read_minim
1924       implicit none
1925       include 'DIMENSIONS'
1926       include 'COMMON.MINIM'
1927       include 'COMMON.IOUNITS'
1928       character*80 ucase
1929       character*320 minimcard
1930       call card_concat(minimcard)
1931       call readi(minimcard,'MAXMIN',maxmin,2000)
1932       call readi(minimcard,'MAXFUN',maxfun,5000)
1933       call readi(minimcard,'MINMIN',minmin,maxmin)
1934       call readi(minimcard,'MINFUN',minfun,maxmin)
1935       call reada(minimcard,'TOLF',tolf,1.0D-2)
1936       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1937       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1938       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1939       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1940       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1941      &         'Options in energy minimization:'
1942       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1943      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1944      & 'MinMin:',MinMin,' MinFun:',MinFun,
1945      & ' TolF:',TolF,' RTolF:',RTolF
1946       return
1947       end
1948 c----------------------------------------------------------------------------
1949       subroutine read_angles(kanal,*)
1950       implicit none
1951       include 'DIMENSIONS'
1952       include 'COMMON.GEO'
1953       include 'COMMON.VAR'
1954       include 'COMMON.CHAIN'
1955       include 'COMMON.IOUNITS'
1956       include 'COMMON.CONTROL'
1957       integer i,kanal
1958 c Read angles from input
1959 c
1960        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1961        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1962        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1963        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1964
1965        do i=1,nres
1966 c 9/7/01 avoid 180 deg valence angle
1967         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1968 c
1969         theta(i)=deg2rad*theta(i)
1970         phi(i)=deg2rad*phi(i)
1971         alph(i)=deg2rad*alph(i)
1972         omeg(i)=deg2rad*omeg(i)
1973        enddo
1974       return
1975    10 return1
1976       end
1977 c----------------------------------------------------------------------------
1978       subroutine reada(rekord,lancuch,wartosc,default)
1979       implicit none
1980       character*(*) rekord,lancuch
1981       double precision wartosc,default
1982       integer ilen,iread
1983       external ilen
1984       iread=index(rekord,lancuch)
1985       if (iread.eq.0) then
1986         wartosc=default
1987         return
1988       endif
1989       iread=iread+ilen(lancuch)+1
1990       read (rekord(iread:),*,err=10,end=10) wartosc
1991       return
1992   10  wartosc=default
1993       return
1994       end
1995 c----------------------------------------------------------------------------
1996       subroutine readi(rekord,lancuch,wartosc,default)
1997       implicit none
1998       character*(*) rekord,lancuch
1999       integer wartosc,default
2000       integer ilen,iread
2001       external ilen
2002       iread=index(rekord,lancuch)
2003       if (iread.eq.0) then
2004         wartosc=default
2005         return
2006       endif
2007       iread=iread+ilen(lancuch)+1
2008       read (rekord(iread:),*,err=10,end=10) wartosc
2009       return
2010   10  wartosc=default
2011       return
2012       end
2013 c----------------------------------------------------------------------------
2014       subroutine multreadi(rekord,lancuch,tablica,dim,default)
2015       implicit none
2016       integer dim,i
2017       integer tablica(dim),default
2018       character*(*) rekord,lancuch
2019       character*80 aux
2020       integer ilen,iread
2021       external ilen
2022       do i=1,dim
2023         tablica(i)=default
2024       enddo
2025       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2026       if (iread.eq.0) return
2027       iread=iread+ilen(lancuch)+1
2028       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2029    10 return
2030       end
2031 c----------------------------------------------------------------------------
2032       subroutine multreada(rekord,lancuch,tablica,dim,default)
2033       implicit none
2034       integer dim,i
2035       double precision tablica(dim),default
2036       character*(*) rekord,lancuch
2037       character*80 aux
2038       integer ilen,iread
2039       external ilen
2040       do i=1,dim
2041         tablica(i)=default
2042       enddo
2043       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2044       if (iread.eq.0) return
2045       iread=iread+ilen(lancuch)+1
2046       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2047    10 return
2048       end
2049 c----------------------------------------------------------------------------
2050       subroutine openunits
2051       implicit none
2052       include 'DIMENSIONS'
2053 #ifdef MPI
2054       include 'mpif.h'
2055       integer ierror
2056       character*16 form,nodename
2057       integer nodelen
2058 #endif
2059       include 'COMMON.SETUP'
2060       include 'COMMON.IOUNITS'
2061       include 'COMMON.MD'
2062       include 'COMMON.CONTROL'
2063       integer lenpre,lenpot,ilen,lentmp,npos
2064       external ilen
2065       character*3 out1file_text,ucase
2066       character*3 ll
2067       external ucase
2068       tmpdir=""
2069       out1file_text="YES"
2070 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2071       call getenv_loc("PREFIX",prefix)
2072       pref_orig = prefix
2073       call getenv_loc("POT",pot)
2074       call getenv_loc("DIRTMP",tmpdir)
2075       call getenv_loc("CURDIR",curdir)
2076       call getenv_loc("OUT1FILE",out1file_text)
2077 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2078       out1file_text=ucase(out1file_text)
2079       if (out1file_text(1:1).eq."Y") then
2080         out1file=.true.
2081       else
2082         out1file=fg_rank.gt.0
2083       endif
2084       lenpre=ilen(prefix)
2085       lenpot=ilen(pot)
2086       lentmp=ilen(tmpdir)
2087       if (lentmp.gt.0) then
2088           write (*,'(80(1h!))')
2089           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2090           write (*,'(80(1h!))')
2091           write (*,*)"All output files will be on node /tmp directory."
2092 #ifdef MPI
2093         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2094         if (me.eq.king) then
2095           write (*,*) "The master node is ",nodename
2096         else if (fg_rank.eq.0) then
2097           write (*,*) "I am the CG slave node ",nodename
2098         else
2099           write (*,*) "I am the FG slave node ",nodename
2100         endif
2101 #endif
2102         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2103         lenpre = lentmp+lenpre+1
2104       endif
2105       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2106 C Get the names and open the input files
2107 #if defined(WINIFL) || defined(WINPGI)
2108       open(1,file=pref_orig(:ilen(pref_orig))//
2109      &  '.inp',status='old',readonly,shared)
2110        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2111 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2112 C Get parameter filenames and open the parameter files.
2113       call getenv_loc('BONDPAR',bondname)
2114       open (ibond,file=bondname,status='old',readonly,shared)
2115       call getenv_loc('THETPAR',thetname)
2116       open (ithep,file=thetname,status='old',readonly,shared)
2117       call getenv_loc('ROTPAR',rotname)
2118       open (irotam,file=rotname,status='old',readonly,shared)
2119       call getenv_loc('TORPAR',torname)
2120       open (itorp,file=torname,status='old',readonly,shared)
2121       call getenv_loc('TORDPAR',tordname)
2122       open (itordp,file=tordname,status='old',readonly,shared)
2123       call getenv_loc('FOURIER',fouriername)
2124       open (ifourier,file=fouriername,status='old',readonly,shared)
2125       call getenv_loc('ELEPAR',elename)
2126       open (ielep,file=elename,status='old',readonly,shared)
2127       call getenv_loc('SIDEPAR',sidename)
2128       open (isidep,file=sidename,status='old',readonly,shared)
2129       call getenv_loc('LIPTRANPAR',liptranname)
2130       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2131 #elif (defined CRAY) || (defined AIX)
2132       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2133      &  action='read')
2134 c      print *,"Processor",myrank," opened file 1"
2135       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2136 c      print *,"Processor",myrank," opened file 9"
2137 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2138 C Get parameter filenames and open the parameter files.
2139       call getenv_loc('BONDPAR',bondname)
2140       open (ibond,file=bondname,status='old',action='read')
2141 c      print *,"Processor",myrank," opened file IBOND"
2142       call getenv_loc('THETPAR',thetname)
2143       open (ithep,file=thetname,status='old',action='read')
2144 #ifndef CRYST_THETA
2145       call getenv_loc('THETPARPDB',thetname_pdb)
2146       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2147 #endif
2148 c      print *,"Processor",myrank," opened file ITHEP"
2149       call getenv_loc('ROTPAR',rotname)
2150       open (irotam,file=rotname,status='old',action='read')
2151 #ifndef CRYST_SC
2152       call getenv_loc('ROTPARPDB',rotname_pdb)
2153       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2154 #endif
2155 c      print *,"Processor",myrank," opened file IROTAM"
2156       call getenv_loc('TORPAR',torname)
2157       open (itorp,file=torname,status='old',action='read')
2158 c      print *,"Processor",myrank," opened file ITORP"
2159       call getenv_loc('TORDPAR',tordname)
2160       open (itordp,file=tordname,status='old',action='read')
2161 c      print *,"Processor",myrank," opened file ITORDP"
2162       call getenv_loc('SCCORPAR',sccorname)
2163       open (isccor,file=sccorname,status='old',action='read')
2164 c      print *,"Processor",myrank," opened file ISCCOR"
2165       call getenv_loc('FOURIER',fouriername)
2166       open (ifourier,file=fouriername,status='old',action='read')
2167 c      print *,"Processor",myrank," opened file IFOURIER"
2168       call getenv_loc('ELEPAR',elename)
2169       open (ielep,file=elename,status='old',action='read')
2170 c      print *,"Processor",myrank," opened file IELEP"
2171       call getenv_loc('SIDEPAR',sidename)
2172       open (isidep,file=sidename,status='old',action='read')
2173       call getenv_loc('LIPTRANPAR',liptranname)
2174       open (iliptranpar,file=liptranname,status='old',action='read')
2175 c      print *,"Processor",myrank," opened file ISIDEP"
2176 c      print *,"Processor",myrank," opened parameter files"
2177 #elif (defined G77)
2178       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2179       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2180 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2181 C Get parameter filenames and open the parameter files.
2182       call getenv_loc('BONDPAR',bondname)
2183       open (ibond,file=bondname,status='old')
2184       call getenv_loc('THETPAR',thetname)
2185       open (ithep,file=thetname,status='old')
2186 #ifndef CRYST_THETA
2187       call getenv_loc('THETPARPDB',thetname_pdb)
2188       open (ithep_pdb,file=thetname_pdb,status='old')
2189 #endif
2190       call getenv_loc('ROTPAR',rotname)
2191       open (irotam,file=rotname,status='old')
2192 #ifndef CRYST_SC
2193       call getenv_loc('ROTPARPDB',rotname_pdb)
2194       open (irotam_pdb,file=rotname_pdb,status='old')
2195 #endif
2196       call getenv_loc('TORPAR',torname)
2197       open (itorp,file=torname,status='old')
2198       call getenv_loc('TORDPAR',tordname)
2199       open (itordp,file=tordname,status='old')
2200       call getenv_loc('SCCORPAR',sccorname)
2201       open (isccor,file=sccorname,status='old')
2202       call getenv_loc('FOURIER',fouriername)
2203       open (ifourier,file=fouriername,status='old')
2204       call getenv_loc('ELEPAR',elename)
2205       open (ielep,file=elename,status='old')
2206       call getenv_loc('SIDEPAR',sidename)
2207       open (isidep,file=sidename,status='old')
2208       call getenv_loc('LIPTRANPAR',liptranname)
2209       open (iliptranpar,file=liptranname,status='old')
2210 #else
2211       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2212      &  readonly)
2213        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2214 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2215 C Get parameter filenames and open the parameter files.
2216       call getenv_loc('BONDPAR',bondname)
2217       open (ibond,file=bondname,status='old',readonly)
2218       call getenv_loc('THETPAR',thetname)
2219       open (ithep,file=thetname,status='old',readonly)
2220       call getenv_loc('ROTPAR',rotname)
2221       open (irotam,file=rotname,status='old',readonly)
2222       call getenv_loc('TORPAR',torname)
2223       open (itorp,file=torname,status='old',readonly)
2224       call getenv_loc('TORDPAR',tordname)
2225       open (itordp,file=tordname,status='old',readonly)
2226       call getenv_loc('SCCORPAR',sccorname)
2227       open (isccor,file=sccorname,status='old',readonly)
2228 #ifndef CRYST_THETA
2229       call getenv_loc('THETPARPDB',thetname_pdb)
2230       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2231 #endif
2232       call getenv_loc('FOURIER',fouriername)
2233       open (ifourier,file=fouriername,status='old',readonly)
2234       call getenv_loc('ELEPAR',elename)
2235       open (ielep,file=elename,status='old',readonly)
2236       call getenv_loc('SIDEPAR',sidename)
2237       open (isidep,file=sidename,status='old',readonly)
2238       call getenv_loc('LIPTRANPAR',liptranname)
2239       open (iliptranpar,file=liptranname,status='old',action='read')
2240 #ifndef CRYST_SC
2241       call getenv_loc('ROTPARPDB',rotname_pdb)
2242       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2243 #endif
2244 #endif
2245 #ifdef TUBE
2246       call getenv_loc('TUBEPAR',tubename)
2247 #if defined(WINIFL) || defined(WINPGI)
2248       open (itube,file=tubename,status='old',readonly,shared)
2249 #elif (defined CRAY)  || (defined AIX)
2250       open (itube,file=tubename,status='old',action='read')
2251 #elif (defined G77)
2252       open (itube,file=tubename,status='old')
2253 #else
2254       open (itube,file=tubename,status='old',readonly)
2255 #endif
2256 #endif
2257 #ifndef OLDSCP
2258 C
2259 C 8/9/01 In the newest version SCp interaction constants are read from a file
2260 C Use -DOLDSCP to use hard-coded constants instead.
2261 C
2262       call getenv_loc('SCPPAR',scpname)
2263 #if defined(WINIFL) || defined(WINPGI)
2264       open (iscpp,file=scpname,status='old',readonly,shared)
2265 #elif (defined CRAY)  || (defined AIX)
2266       open (iscpp,file=scpname,status='old',action='read')
2267 #elif (defined G77)
2268       open (iscpp,file=scpname,status='old')
2269 #else
2270       open (iscpp,file=scpname,status='old',readonly)
2271 #endif
2272 #endif
2273       call getenv_loc('PATTERN',patname)
2274 #if defined(WINIFL) || defined(WINPGI)
2275       open (icbase,file=patname,status='old',readonly,shared)
2276 #elif (defined CRAY)  || (defined AIX)
2277       open (icbase,file=patname,status='old',action='read')
2278 #elif (defined G77)
2279       open (icbase,file=patname,status='old')
2280 #else
2281       open (icbase,file=patname,status='old',readonly)
2282 #endif
2283 #ifdef MPI
2284 C Open output file only for CG processes
2285 c      print *,"Processor",myrank," fg_rank",fg_rank
2286       if (fg_rank.eq.0) then
2287
2288       if (nodes.eq.1) then
2289         npos=3
2290       else
2291         npos = dlog10(dfloat(nodes-1))+1
2292       endif
2293       if (npos.lt.3) npos=3
2294       write (liczba,'(i1)') npos
2295       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2296      &  //')'
2297       write (liczba,form) me
2298       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2299      &  liczba(:ilen(liczba))
2300       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2301      &  //'.int'
2302       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2303      &  //'.pdb'
2304       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2305      &  liczba(:ilen(liczba))//'.mol2'
2306       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2307      &  liczba(:ilen(liczba))//'.stat'
2308       if (lentmp.gt.0)
2309      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2310      &      //liczba(:ilen(liczba))//'.stat')
2311       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2312      &  //'.rst'
2313       if(usampl) then
2314           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2315      & liczba(:ilen(liczba))//'.const'
2316       endif
2317
2318       endif
2319 #else
2320       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2321       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2322       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2323       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2324       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2325       if (lentmp.gt.0)
2326      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2327      &    //'.stat')
2328       rest2name=prefix(:ilen(prefix))//'.rst'
2329       if(usampl) then
2330          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2331       endif
2332 #endif
2333 #if defined(AIX) || defined(PGI) || defined(CRAY)
2334       if (me.eq.king .or. .not. out1file) then
2335          open(iout,file=outname,status='unknown')
2336       else
2337          open(iout,file="/dev/null",status="unknown")
2338       endif
2339 c#define DEBUG
2340 #ifdef DEBUG
2341       if (fg_rank.gt.0) then
2342         write (liczba,'(i3.3)') myrank/nfgtasks
2343         write (ll,'(bz,i3.3)') fg_rank
2344         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2345      &   status='unknown')
2346       endif
2347 #endif
2348 c#undef DEBUG
2349       if(me.eq.king) then
2350        open(igeom,file=intname,status='unknown',position='append')
2351        open(ipdb,file=pdbname,status='unknown')
2352        open(imol2,file=mol2name,status='unknown')
2353        open(istat,file=statname,status='unknown',position='append')
2354       else
2355 c1out       open(iout,file=outname,status='unknown')
2356       endif
2357 #else
2358       if (me.eq.king .or. .not.out1file)
2359      &    open(iout,file=outname,status='unknown')
2360 #define DEBUG
2361 #ifdef DEBUG
2362       if (fg_rank.gt.0) then
2363         write (liczba,'(i3.3)') myrank/nfgtasks
2364         write (ll,'(bz,i3.3)') fg_rank
2365         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2366      &   status='unknown')
2367       endif
2368 #endif
2369 #undef DEBUG
2370       if(me.eq.king) then
2371        open(igeom,file=intname,status='unknown',access='append')
2372        open(ipdb,file=pdbname,status='unknown')
2373        open(imol2,file=mol2name,status='unknown')
2374        open(istat,file=statname,status='unknown',access='append')
2375       else
2376 c1out       open(iout,file=outname,status='unknown')
2377       endif
2378 #endif
2379       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2380       csa_seed=prefix(:lenpre)//'.CSA.seed'
2381       csa_history=prefix(:lenpre)//'.CSA.history'
2382       csa_bank=prefix(:lenpre)//'.CSA.bank'
2383       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2384       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2385       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2386 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2387       csa_int=prefix(:lenpre)//'.int'
2388       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2389       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2390       csa_in=prefix(:lenpre)//'.CSA.in'
2391 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2392 C Write file names
2393       if (me.eq.king)then
2394       write (iout,'(80(1h-))')
2395       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2396       write (iout,'(80(1h-))')
2397       write (iout,*) "Input file                      : ",
2398      &  pref_orig(:ilen(pref_orig))//'.inp'
2399       write (iout,*) "Output file                     : ",
2400      &  outname(:ilen(outname))
2401       write (iout,*)
2402       write (iout,*) "Sidechain potential file        : ",
2403      &  sidename(:ilen(sidename))
2404 #ifndef OLDSCP
2405       write (iout,*) "SCp potential file              : ",
2406      &  scpname(:ilen(scpname))
2407 #endif
2408       write (iout,*) "Electrostatic potential file    : ",
2409      &  elename(:ilen(elename))
2410       write (iout,*) "Cumulant coefficient file       : ",
2411      &  fouriername(:ilen(fouriername))
2412       write (iout,*) "Torsional parameter file        : ",
2413      &  torname(:ilen(torname))
2414       write (iout,*) "Double torsional parameter file : ",
2415      &  tordname(:ilen(tordname))
2416       write (iout,*) "SCCOR parameter file : ",
2417      &  sccorname(:ilen(sccorname))
2418       write (iout,*) "Bond & inertia constant file    : ",
2419      &  bondname(:ilen(bondname))
2420       write (iout,*) "Bending-torsion parameter file  : ",
2421      &  thetname(:ilen(thetname))
2422       write (iout,*) "Rotamer parameter file          : ",
2423      &  rotname(:ilen(rotname))
2424       write (iout,*) "Thetpdb parameter file          : ",
2425      &  thetname_pdb(:ilen(thetname_pdb))
2426       write (iout,*) "Threading database              : ",
2427      &  patname(:ilen(patname))
2428       if (lentmp.ne.0)
2429      &write (iout,*)" DIRTMP                          : ",
2430      &  tmpdir(:lentmp)
2431       write (iout,'(80(1h-))')
2432       endif
2433       return
2434       end
2435 c----------------------------------------------------------------------------
2436       subroutine card_concat(card)
2437       implicit real*8 (a-h,o-z)
2438       include 'DIMENSIONS'
2439       include 'COMMON.IOUNITS'
2440       character*(*) card
2441       character*80 karta,ucase
2442       external ilen
2443       read (inp,'(a)') karta
2444       karta=ucase(karta)
2445       card=' '
2446       do while (karta(80:80).eq.'&')
2447         card=card(:ilen(card)+1)//karta(:79)
2448         read (inp,'(a)') karta
2449         karta=ucase(karta)
2450       enddo
2451       card=card(:ilen(card)+1)//karta
2452       return
2453       end
2454 c------------------------------------------------------------------------------
2455       subroutine readrst
2456       implicit none
2457       include 'DIMENSIONS'
2458       include 'COMMON.CHAIN'
2459       include 'COMMON.IOUNITS'
2460       include 'COMMON.CONTROL'
2461       include 'COMMON.MD'
2462       include 'COMMON.QRESTR'
2463 #ifdef FIVEDIAG
2464        include 'COMMON.LAGRANGE.5diag'
2465 #else
2466        include 'COMMON.LAGRANGE'
2467 #endif
2468       integer i,j
2469       open(irest2,file=rest2name,status='unknown')
2470       read(irest2,*) totT,EK,potE,totE,t_bath
2471       totTafm=totT
2472       do i=0,2*nres-1
2473          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2474       enddo
2475       do i=0,2*nres-1
2476          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2477       enddo
2478       if(usampl) then
2479              read (irest2,*) iset
2480       endif
2481       close(irest2)
2482       return
2483       end
2484 c------------------------------------------------------------------------------
2485       subroutine read_fragments
2486       implicit none
2487       include 'DIMENSIONS'
2488 #ifdef MPI
2489       include 'mpif.h'
2490 #endif
2491       include 'COMMON.SETUP'
2492       include 'COMMON.CHAIN'
2493       include 'COMMON.IOUNITS'
2494       include 'COMMON.MD'
2495       include 'COMMON.QRESTR'
2496       include 'COMMON.CONTROL'
2497       integer i
2498       read(inp,*) nset,nfrag,npair,nfrag_back
2499       loc_qlike=(nfrag_back.lt.0)
2500       nfrag_back=iabs(nfrag_back)
2501       if(me.eq.king.or..not.out1file)
2502      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2503      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2504       do iset=1,nset
2505          read(inp,*) mset(iset)
2506        do i=1,nfrag
2507          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2508      &     qinfrag(i,iset)
2509          if(me.eq.king.or..not.out1file)
2510      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2511      &     ifrag(2,i,iset), qinfrag(i,iset)
2512        enddo
2513        do i=1,npair
2514         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2515      &    qinpair(i,iset)
2516         if(me.eq.king.or..not.out1file)
2517      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2518      &    ipair(2,i,iset), qinpair(i,iset)
2519        enddo
2520        if (loc_qlike) then
2521        do i=1,nfrag_back
2522         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2523      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2524      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2525      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2526         if(me.eq.king.or..not.out1file)
2527      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2528      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2529      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2530      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2531        enddo
2532        else
2533        do i=1,nfrag_back
2534         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2535      &     wfrag_back(3,i,iset),
2536      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2537         if(me.eq.king.or..not.out1file)
2538      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2539      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2540        enddo
2541        endif
2542       enddo
2543       return
2544       end
2545 C---------------------------------------------------------------------------
2546       subroutine read_afminp
2547       implicit none
2548       include 'DIMENSIONS'
2549 #ifdef MPI
2550       include 'mpif.h'
2551 #endif
2552       include 'COMMON.SETUP'
2553       include 'COMMON.CONTROL'
2554       include 'COMMON.CHAIN'
2555       include 'COMMON.IOUNITS'
2556       include 'COMMON.SBRIDGE'
2557       character*320 afmcard
2558       integer i
2559 c      print *, "wchodze"
2560       call card_concat(afmcard)
2561       call readi(afmcard,"BEG",afmbeg,0)
2562       call readi(afmcard,"END",afmend,0)
2563       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2564       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2565 c      print *,'FORCE=' ,forceAFMconst
2566 CCCC NOW PROPERTIES FOR AFM
2567        distafminit=0.0d0
2568        do i=1,3
2569         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2570        enddo
2571         distafminit=dsqrt(distafminit)
2572 c        print *,'initdist',distafminit
2573       return
2574       end
2575 c-------------------------------------------------------------------------------
2576       subroutine read_saxs_constr
2577       implicit none
2578       include 'DIMENSIONS'
2579 #ifdef MPI
2580       include 'mpif.h'
2581 #endif
2582       include 'COMMON.SETUP'
2583       include 'COMMON.CONTROL'
2584       include 'COMMON.SAXS'
2585       include 'COMMON.CHAIN'
2586       include 'COMMON.IOUNITS'
2587       include 'COMMON.SBRIDGE'
2588       double precision cm(3),cnorm
2589       integer i,j
2590 c      read(inp,*) nsaxs
2591       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2592       call flush(iout)
2593       if (saxs_mode.eq.0) then
2594 c SAXS distance distribution
2595       do i=1,nsaxs
2596         read(inp,*) distsaxs(i),Psaxs(i)
2597       enddo
2598       Cnorm = 0.0d0
2599       do i=1,nsaxs
2600         Cnorm = Cnorm + Psaxs(i)
2601       enddo
2602       write (iout,*) "Cnorm",Cnorm
2603       do i=1,nsaxs
2604         Psaxs(i)=Psaxs(i)/Cnorm
2605       enddo
2606       write (iout,*) "Normalized distance distribution from SAXS"
2607       do i=1,nsaxs
2608         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2609       enddo
2610       Wsaxs0=0.0d0
2611       do i=1,nsaxs
2612         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2613       enddo
2614       write (iout,*) "Wsaxs0",Wsaxs0
2615       else
2616 c SAXS "spheres".
2617       do i=1,nsaxs
2618         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2619       enddo
2620       do j=1,3
2621         cm(j)=0.0d0
2622       enddo
2623       do i=1,nsaxs
2624         do j=1,3
2625           cm(j)=cm(j)+Csaxs(j,i)
2626         enddo
2627       enddo
2628       do j=1,3
2629         cm(j)=cm(j)/nsaxs
2630       enddo
2631       do i=1,nsaxs
2632         do j=1,3
2633           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2634         enddo
2635       enddo
2636       write (iout,*) "SAXS sphere coordinates"
2637       do i=1,nsaxs
2638         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2639       enddo
2640       endif
2641       return
2642       end
2643
2644 c-------------------------------------------------------------------------------
2645       subroutine read_dist_constr
2646       implicit none
2647       include 'DIMENSIONS'
2648 #ifdef MPI
2649       include 'mpif.h'
2650 #endif
2651       include 'COMMON.SETUP'
2652       include 'COMMON.CONTROL'
2653       include 'COMMON.CHAIN'
2654       include 'COMMON.IOUNITS'
2655       include 'COMMON.SBRIDGE'
2656       include 'COMMON.INTERACT'
2657       integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2658       integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2659       double precision wfrag_(100),wpair_(1000)
2660       double precision ddjk,dist,dist_cut,fordepthmax
2661       character*5000 controlcard
2662       logical normalize,next
2663       integer restr_type
2664       double precision scal_bfac
2665       double precision xlink(4,0:4) /
2666 c           a          b       c     sigma
2667      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2668      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2669      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2670      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2671      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2672 c      print *, "WCHODZE"
2673       write (iout,*) "Calling read_dist_constr"
2674 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2675 c      call flush(iout)
2676       restr_on_coord=.false.
2677       next=.true.
2678
2679       npeak=0
2680       ipeak=0
2681       nhpb_peak=0
2682
2683       DO WHILE (next)
2684
2685       call card_concat(controlcard)
2686       next = index(controlcard,"NEXT").gt.0
2687       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2688       write (iout,*) "restr_type",restr_type
2689       call readi(controlcard,"NFRAG",nfrag_,0)
2690       call readi(controlcard,"NFRAG",nfrag_,0)
2691       call readi(controlcard,"NPAIR",npair_,0)
2692       call readi(controlcard,"NDIST",ndist_,0)
2693       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2694       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2695       if (restr_type.eq.10)
2696      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2697       if (restr_type.eq.12)
2698      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2699       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2700       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2701       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2702       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2703       normalize = index(controlcard,"NORMALIZE").gt.0
2704       write (iout,*) "WBOLTZD",wboltzd
2705       write (iout,*) "SCAL_PEAK",scal_peak
2706       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2707       write (iout,*) "IFRAG"
2708       do i=1,nfrag_
2709         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2710       enddo
2711       write (iout,*) "IPAIR"
2712       do i=1,npair_
2713         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2714       enddo
2715       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2716      & write (iout,*)
2717      &   "Distance restraints as generated from reference structure"
2718       do i=1,nfrag_
2719         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2720         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2721      &    ifrag_(2,i)=nstart_sup+nsup-1
2722 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2723 c        call flush(iout)
2724         if (wfrag_(i).eq.0.0d0) cycle
2725         do j=ifrag_(1,i),ifrag_(2,i)-1
2726           do k=j+1,ifrag_(2,i)
2727 c            write (iout,*) "j",j," k",k
2728             ddjk=dist(j,k)
2729             if (restr_type.eq.1) then
2730               nhpb=nhpb+1
2731               irestr_type(nhpb)=1
2732               ihpb(nhpb)=j
2733               jhpb(nhpb)=k
2734               dhpb(nhpb)=ddjk
2735               forcon(nhpb)=wfrag_(i)
2736             else if (constr_dist.eq.2) then
2737               if (ddjk.le.dist_cut) then
2738                 nhpb=nhpb+1
2739                 irestr_type(nhpb)=1
2740                 ihpb(nhpb)=j
2741                 jhpb(nhpb)=k
2742                 dhpb(nhpb)=ddjk
2743                 forcon(nhpb)=wfrag_(i)
2744               endif
2745             else if (restr_type.eq.3) then
2746               nhpb=nhpb+1
2747               irestr_type(nhpb)=1
2748               ihpb(nhpb)=j
2749               jhpb(nhpb)=k
2750               dhpb(nhpb)=ddjk
2751               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2752             endif
2753 #ifdef MPI
2754             if (.not.out1file .or. me.eq.king)
2755      &      write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
2756      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2757 #else
2758             write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
2759      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2760 #endif
2761           enddo
2762         enddo
2763       enddo
2764       do i=1,npair_
2765         if (wpair_(i).eq.0.0d0) cycle
2766         ii = ipair_(1,i)
2767         jj = ipair_(2,i)
2768         if (ii.gt.jj) then
2769           itemp=ii
2770           ii=jj
2771           jj=itemp
2772         endif
2773         do j=ifrag_(1,ii),ifrag_(2,ii)
2774           do k=ifrag_(1,jj),ifrag_(2,jj)
2775             ddjk=dist(j,k)
2776             if (restr_type.eq.1) then
2777               nhpb=nhpb+1
2778               irestr_type(nhpb)=1
2779               ihpb(nhpb)=j
2780               jhpb(nhpb)=k
2781               dhpb(nhpb)=ddjk
2782               forcon(nhpb)=wpair_(i)
2783             else if (constr_dist.eq.2) then
2784               if (ddjk.le.dist_cut) then
2785                 nhpb=nhpb+1
2786                 irestr_type(nhpb)=1
2787                 ihpb(nhpb)=j
2788                 jhpb(nhpb)=k
2789                 dhpb(nhpb)=ddjk
2790                 forcon(nhpb)=wpair_(i)
2791               endif
2792             else if (restr_type.eq.3) then
2793               nhpb=nhpb+1
2794               irestr_type(nhpb)=1
2795               ihpb(nhpb)=j
2796               jhpb(nhpb)=k
2797               dhpb(nhpb)=ddjk
2798               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2799             endif
2800 #ifdef MPI
2801             if (.not.out1file .or. me.eq.king)
2802      &      write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
2803      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2804 #else
2805             write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
2806      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2807 #endif
2808           enddo
2809         enddo
2810       enddo
2811
2812 c      print *,ndist_
2813       write (iout,*) "Distance restraints as read from input"
2814       do i=1,ndist_
2815         if (restr_type.eq.12) then
2816           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2817      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2818      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2819      &    fordepth_peak(nhpb_peak+1),npeak
2820 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2821 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2822 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2823 c     &    fordepth_peak(nhpb_peak+1),npeak
2824           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2825      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2826           nhpb_peak=nhpb_peak+1
2827           irestr_type_peak(nhpb_peak)=12
2828           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2829           ipeak(2,npeak)=i
2830 #ifdef MPI
2831           if (.not.out1file .or. me.eq.king)
2832      &    write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
2833      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2834      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2835      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2836      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2837 #else
2838           write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
2839      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2840      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2841      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2842      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2843 #endif
2844           if (ibecarb_peak(nhpb_peak).eq.3) then
2845             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2846           else if (ibecarb_peak(nhpb_peak).eq.2) then
2847             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2848           else if (ibecarb_peak(nhpb_peak).eq.1) then
2849             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2850             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2851           endif
2852         else if (restr_type.eq.11) then
2853           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2854      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2855 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2856           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2857           nhpb=nhpb+1
2858           irestr_type(nhpb)=11
2859 #ifdef MPI
2860           if (.not.out1file .or. me.eq.king)
2861      &    write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
2862      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2863      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2864 #else
2865           write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
2866      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2867      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2868 #endif
2869 c          if (ibecarb(nhpb).gt.0) then
2870 c            ihpb(nhpb)=ihpb(nhpb)+nres
2871 c            jhpb(nhpb)=jhpb(nhpb)+nres
2872 c          endif
2873          if (ibecarb(nhpb).eq.3) then
2874             ihpb(nhpb)=ihpb(nhpb)+nres
2875           else if (ibecarb(nhpb).eq.2) then
2876             ihpb(nhpb)=ihpb(nhpb)+nres
2877           else if (ibecarb(nhpb).eq.1) then
2878             ihpb(nhpb)=ihpb(nhpb)+nres
2879             jhpb(nhpb)=jhpb(nhpb)+nres
2880           endif
2881         else if (restr_type.eq.10) then
2882 c Cross-lonk Markov-like potential
2883           call card_concat(controlcard)
2884           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2885           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2886           ibecarb(nhpb+1)=0
2887           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2888           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2889           if (index(controlcard,"ZL").gt.0) then
2890             link_type=1
2891           else if (index(controlcard,"ADH").gt.0) then
2892             link_type=2
2893           else if (index(controlcard,"PDH").gt.0) then
2894             link_type=3
2895           else if (index(controlcard,"DSS").gt.0) then
2896             link_type=4
2897           else
2898             link_type=0
2899           endif
2900           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2901      &       xlink(1,link_type))
2902           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2903      &       xlink(2,link_type))
2904           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2905      &       xlink(3,link_type))
2906           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2907      &       xlink(4,link_type))
2908           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2909 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2910 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2911           if (forcon(nhpb+1).le.0.0d0 .or.
2912      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2913           nhpb=nhpb+1
2914           irestr_type(nhpb)=10
2915           if (ibecarb(nhpb).eq.3) then
2916             jhpb(nhpb)=jhpb(nhpb)+nres
2917           else if (ibecarb(nhpb).eq.2) then
2918             ihpb(nhpb)=ihpb(nhpb)+nres
2919           else if (ibecarb(nhpb).eq.1) then
2920             ihpb(nhpb)=ihpb(nhpb)+nres
2921             jhpb(nhpb)=jhpb(nhpb)+nres
2922           endif
2923 #ifdef MPI
2924           if (.not.out1file .or. me.eq.king)
2925      &    write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
2926      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2927      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2928      &     irestr_type(nhpb)
2929 #else
2930           write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
2931      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2932      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2933      &     irestr_type(nhpb)
2934 #endif
2935         else
2936 C        print *,"in else"
2937           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2938      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2939           if (forcon(nhpb+1).gt.0.0d0) then
2940           nhpb=nhpb+1
2941           if (dhpb1(nhpb).eq.0.0d0) then
2942             irestr_type(nhpb)=1
2943           else
2944             irestr_type(nhpb)=2
2945           endif
2946           if (ibecarb(nhpb).eq.3) then
2947             jhpb(nhpb)=jhpb(nhpb)+nres
2948           else if (ibecarb(nhpb).eq.2) then
2949             ihpb(nhpb)=ihpb(nhpb)+nres
2950           else if (ibecarb(nhpb).eq.1) then
2951             ihpb(nhpb)=ihpb(nhpb)+nres
2952             jhpb(nhpb)=jhpb(nhpb)+nres
2953           endif
2954           if (dhpb(nhpb).eq.0.0d0)
2955      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2956         endif
2957 #ifdef MPI
2958           if (.not.out1file .or. me.eq.king)
2959      &    write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
2960      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2961 #else
2962           write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
2963      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2964 #endif
2965         endif
2966 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2967 C        if (forcon(nhpb+1).gt.0.0d0) then
2968 C          nhpb=nhpb+1
2969 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2970       enddo
2971
2972       if (restr_type.eq.4) then
2973         write (iout,*) "The BFAC array"
2974         do i=nnt,nct
2975           write (iout,'(i5,f10.5)') i,bfac(i)
2976         enddo
2977         do i=nnt,nct
2978           if (itype(i).eq.ntyp1) cycle
2979           do j=nnt,i-1
2980             if (itype(j).eq.ntyp1) cycle
2981             if (itype(i).eq.10) then
2982               iiend=0
2983             else
2984               iiend=1
2985             endif
2986             if (itype(j).eq.10) then
2987               jjend=0
2988             else
2989               jjend=1
2990             endif
2991             kk=0
2992             do ii=0,iiend
2993             do jj=0,jjend
2994             nhpb=nhpb+1
2995             irestr_type(nhpb)=1
2996             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2997             irestr_type(nhpb)=1
2998             ibecarb(nhpb)=kk
2999             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
3000             ihpb(nhpb)=i+nres*ii
3001             jhpb(nhpb)=j+nres*jj
3002             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
3003 #ifdef MPI
3004             if (.not.out1file .or. me.eq.king) then
3005             write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
3006      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
3007      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
3008      &       irestr_type(nhpb)
3009             endif
3010 #else
3011             write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
3012      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
3013      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
3014      &       irestr_type(nhpb)
3015 #endif
3016             kk=kk+1
3017           enddo
3018           enddo
3019           enddo
3020         enddo
3021       endif
3022
3023       if (restr_type.eq.5) then
3024         restr_on_coord=.true.
3025         do i=nnt,nct
3026           if (itype(i).eq.ntyp1) cycle
3027           bfac(i)=(scal_bfac/bfac(i))**2
3028         enddo
3029       endif
3030
3031       ENDDO ! next
3032
3033       fordepthmax=0.0d0
3034       if (normalize) then
3035         do i=nss+1,nhpb
3036           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3037      &      fordepthmax=fordepth(i)
3038         enddo
3039         do i=nss+1,nhpb
3040           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3041         enddo
3042       endif
3043       return
3044       end
3045 c-------------------------------------------------------------------------------
3046       subroutine read_constr_homology
3047       implicit none
3048       include 'DIMENSIONS'
3049 #ifdef MPI
3050       include 'mpif.h'
3051 #endif
3052       include 'COMMON.SETUP'
3053       include 'COMMON.CONTROL'
3054       include 'COMMON.HOMOLOGY'
3055       include 'COMMON.CHAIN'
3056       include 'COMMON.IOUNITS'
3057       include 'COMMON.MD'
3058       include 'COMMON.QRESTR'
3059       include 'COMMON.GEO'
3060       include 'COMMON.INTERACT'
3061       include 'COMMON.NAMES'
3062 c
3063 c For new homol impl
3064 c
3065       include 'COMMON.VAR'
3066 c
3067
3068 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3069 c    &                 dist_cut
3070 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3071 c    &    sigma_odl_temp(maxres,maxres,max_template)
3072       character*2 kic2
3073       character*24 model_ki_dist, model_ki_angle
3074       character*500 controlcard
3075       integer ki,i,ii,j,k,l,ii_in_use(maxchain*maxdim),i_tmp,
3076      & idomain_tmp,
3077      & irec,ik,iistart,nres_temp
3078       integer ilen
3079       external ilen
3080       logical liiflag,lfirst
3081       integer i01,i10
3082 c
3083 c     FP - Nov. 2014 Temporary specifications for new vars
3084 c
3085       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3086      &                 rescore3_tmp
3087       double precision, dimension (max_template,maxres) :: rescore
3088       double precision, dimension (max_template,maxres) :: rescore2
3089       double precision, dimension (max_template,maxres) :: rescore3
3090       double precision distal
3091       character*24 tpl_k_rescore
3092       character*256 pdbfile
3093 c -----------------------------------------------------------------
3094 c Reading multiple PDB ref structures and calculation of retraints
3095 c not using pre-computed ones stored in files model_ki_{dist,angle}
3096 c FP (Nov., 2014)
3097 c -----------------------------------------------------------------
3098 c
3099 c
3100 c Alternative: reading from input
3101       call card_concat(controlcard)
3102       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3103       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3104       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3105       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3106       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3107       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3108       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3109       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3110       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3111       if(.not.read2sigma.and.start_from_model) then
3112           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3113      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3114           start_from_model=.false.
3115           iranconf=(indpdb.le.0)
3116       endif
3117       if(start_from_model .and. (me.eq.king .or. .not. out1file))
3118      &    write(iout,*) 'START_FROM_MODELS is ON'
3119 c      if(start_from_model .and. rest) then
3120 c        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3121 c         write(iout,*) 'START_FROM_MODELS is OFF'
3122 c         write(iout,*) 'remove restart keyword from input'
3123 c        endif
3124 c      endif
3125       if (start_from_model) nmodel_start=constr_homology
3126       if (homol_nset.gt.1)then
3127          call card_concat(controlcard)
3128          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3129          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3130 c          write(iout,*) "iset homology_weight "
3131           do i=1,homol_nset
3132            write(iout,*) i,waga_homology(i)
3133           enddo
3134          endif
3135          iset=mod(kolor,homol_nset)+1
3136       else
3137        iset=1
3138        waga_homology(1)=1.0
3139       endif
3140
3141 cd      write (iout,*) "nnt",nnt," nct",nct
3142 cd      call flush(iout)
3143
3144       if (read_homol_frag) then
3145         call read_klapaucjusz
3146       else
3147
3148       lim_odl=0
3149       lim_dih=0
3150 c
3151 c      write(iout,*) 'nnt=',nnt,'nct=',nct
3152 c
3153 c      do i = nnt,nct
3154 c        do k=1,constr_homology
3155 c         idomain(k,i)=0
3156 c        enddo
3157 c      enddo
3158        idomain=0
3159
3160 c      ii=0
3161 c      do i = nnt,nct-2
3162 c        do j=i+2,nct
3163 c        ii=ii+1
3164 c        ii_in_use(ii)=0
3165 c        enddo
3166 c      enddo
3167       ii_in_use=0
3168       do k=1,constr_homology
3169
3170         read(inp,'(a)') pdbfile
3171         pdbfiles_chomo(k)=pdbfile
3172         if(me.eq.king .or. .not. out1file)
3173      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3174      &  pdbfile(:ilen(pdbfile))
3175         open(ipdbin,file=pdbfile,status='old',err=33)
3176         goto 34
3177   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3178      &  pdbfile(:ilen(pdbfile))
3179         stop
3180   34    continue
3181 c        print *,'Begin reading pdb data'
3182 c
3183 c Files containing res sim or local scores (former containing sigmas)
3184 c
3185
3186         write(kic2,'(bz,i2.2)') k
3187
3188         tpl_k_rescore="template"//kic2//".sco"
3189
3190         unres_pdb=.false.
3191         nres_temp=nres
3192         if (read2sigma) then
3193           call readpdb_template(k)
3194         else
3195           call readpdb
3196         endif
3197         nres_chomo(k)=nres
3198         nres=nres_temp
3199 c
3200 c     Distance restraints
3201 c
3202 c          ... --> odl(k,ii)
3203 C Copy the coordinates from reference coordinates (?)
3204         do i=1,2*nres_chomo(k)
3205           do j=1,3
3206             c(j,i)=cref(j,i)
3207 c           write (iout,*) "c(",j,i,") =",c(j,i)
3208           enddo
3209         enddo
3210 c
3211 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3212 c
3213 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3214           open (ientin,file=tpl_k_rescore,status='old')
3215           if (nnt.gt.1) rescore(k,1)=0.0d0
3216           do irec=nnt,nct ! loop for reading res sim
3217             if (read2sigma) then
3218              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3219      &                                rescore3_tmp,idomain_tmp
3220              i_tmp=i_tmp+nnt-1
3221              idomain(k,i_tmp)=idomain_tmp
3222              rescore(k,i_tmp)=rescore_tmp
3223              rescore2(k,i_tmp)=rescore2_tmp
3224              rescore3(k,i_tmp)=rescore3_tmp
3225              if (.not. out1file .or. me.eq.king)
3226      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3227      &                      i_tmp,rescore2_tmp,rescore_tmp,
3228      &                                rescore3_tmp,idomain_tmp
3229             else
3230              idomain(k,irec)=1
3231              read (ientin,*,end=1401) rescore_tmp
3232
3233 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3234              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3235 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3236             endif
3237           enddo
3238  1401   continue
3239         close (ientin)
3240         if (waga_dist.ne.0.0d0) then
3241           ii=0
3242           do i = nnt,nct-2
3243             do j=i+2,nct
3244
3245               x12=c(1,i)-c(1,j)
3246               y12=c(2,i)-c(2,j)
3247               z12=c(3,i)-c(3,j)
3248               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3249 c              write (iout,*) k,i,j,distal,dist2_cut
3250
3251             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3252      &            .and. distal.le.dist2_cut ) then
3253
3254               ii=ii+1
3255               ii_in_use(ii)=1
3256               l_homo(k,ii)=.true.
3257
3258 c             write (iout,*) "k",k
3259 c             write (iout,*) "i",i," j",j," constr_homology",
3260 c    &                       constr_homology
3261               ires_homo(ii)=i
3262               jres_homo(ii)=j
3263               odl(k,ii)=distal
3264               if (read2sigma) then
3265                 sigma_odl(k,ii)=0
3266                 do ik=i,j
3267                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3268                 enddo
3269                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3270                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3271      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3272               else
3273                 if (odl(k,ii).le.dist_cut) then
3274                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3275                 else
3276 #ifdef OLDSIGMA
3277                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3278      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3279 #else
3280                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3281      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3282 #endif
3283                 endif
3284               endif
3285               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3286             else
3287 c              ii=ii+1
3288 c              l_homo(k,ii)=.false.
3289             endif
3290             enddo
3291           enddo
3292         lim_odl=ii
3293         endif
3294 c        write (iout,*) "Distance restraints set"
3295 c        call flush(iout)
3296 c
3297 c     Theta, dihedral and SC retraints
3298 c
3299         if (waga_angle.gt.0.0d0) then
3300 c         open (ientin,file=tpl_k_sigma_dih,status='old')
3301 c         do irec=1,maxres-3 ! loop for reading sigma_dih
3302 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3303 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3304 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3305 c    &                            sigma_dih(k,i+nnt-1)
3306 c         enddo
3307 c1402   continue
3308 c         close (ientin)
3309           do i = nnt+3,nct
3310             if (idomain(k,i).eq.0) then
3311                sigma_dih(k,i)=0.0
3312                cycle
3313             endif
3314             dih(k,i)=phiref(i) ! right?
3315 c           read (ientin,*) sigma_dih(k,i) ! original variant
3316 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
3317 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3318 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3319 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
3320 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
3321
3322             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3323      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3324 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3325 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
3326 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3327 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3328 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
3329 c   Instead of res sim other local measure of b/b str reliability possible
3330             if (sigma_dih(k,i).ne.0)
3331      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3332 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3333           enddo
3334           lim_dih=nct-nnt-2
3335         endif
3336 c        write (iout,*) "Dihedral angle restraints set"
3337 c        call flush(iout)
3338
3339         if (waga_theta.gt.0.0d0) then
3340 c         open (ientin,file=tpl_k_sigma_theta,status='old')
3341 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3342 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3343 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3344 c    &                              sigma_theta(k,i+nnt-1)
3345 c         enddo
3346 c1403   continue
3347 c         close (ientin)
3348
3349           do i = nnt+2,nct ! right? without parallel.
3350 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
3351 c         do i=ithet_start,ithet_end ! with FG parallel.
3352              if (idomain(k,i).eq.0) then
3353               sigma_theta(k,i)=0.0
3354               cycle
3355              endif
3356              thetatpl(k,i)=thetaref(i)
3357 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3358 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
3359 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3360 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
3361 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
3362              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3363      &                        rescore(k,i-2))/3.0
3364 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3365              if (sigma_theta(k,i).ne.0)
3366      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3367
3368 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3369 c                             rescore(k,i-2) !  right expression ?
3370 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3371           enddo
3372         endif
3373 c        write (iout,*) "Angle restraints set"
3374 c        call flush(iout)
3375
3376         if (waga_d.gt.0.0d0) then
3377 c       open (ientin,file=tpl_k_sigma_d,status='old')
3378 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3379 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3380 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3381 c    &                          sigma_d(k,i+nnt-1)
3382 c         enddo
3383 c1404   continue
3384
3385           do i = nnt,nct ! right? without parallel.
3386 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3387 c         do i=loc_start,loc_end ! with FG parallel.
3388                if (itype(i).eq.10) cycle
3389                if (idomain(k,i).eq.0 ) then
3390                   sigma_d(k,i)=0.0
3391                   cycle
3392                endif
3393                xxtpl(k,i)=xxref(i)
3394                yytpl(k,i)=yyref(i)
3395                zztpl(k,i)=zzref(i)
3396 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3397 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3398 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3399 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3400                sigma_d(k,i)=rescore3(k,i) !  right expression ?
3401                if (sigma_d(k,i).ne.0)
3402      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3403
3404 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3405 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3406 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3407           enddo
3408         endif
3409       enddo
3410 c      write (iout,*) "SC restraints set"
3411 c      call flush(iout)
3412 c
3413 c remove distance restraints not used in any model from the list
3414 c shift data in all arrays
3415 c
3416 c      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3417       if (waga_dist.ne.0.0d0) then
3418         ii=0
3419         liiflag=.true.
3420         lfirst=.true.
3421         do i=nnt,nct-2
3422          do j=i+2,nct
3423           ii=ii+1
3424 c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3425 c     &            .and. distal.le.dist2_cut ) then
3426 c          write (iout,*) "i",i," j",j," ii",ii
3427 c          call flush(iout)
3428           if (ii_in_use(ii).eq.0.and.liiflag.or.
3429      &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3430             liiflag=.false.
3431             i10=ii
3432             if (lfirst) then
3433               lfirst=.false.
3434               iistart=ii
3435             else
3436               if(i10.eq.lim_odl) i10=i10+1
3437               do ki=0,i10-i01-1
3438                ires_homo(iistart+ki)=ires_homo(ki+i01)
3439                jres_homo(iistart+ki)=jres_homo(ki+i01)
3440                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3441                do k=1,constr_homology
3442                 odl(k,iistart+ki)=odl(k,ki+i01)
3443                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3444                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3445                enddo
3446               enddo
3447               iistart=iistart+i10-i01
3448             endif
3449           endif
3450           if (ii_in_use(ii).ne.0.and..not.liiflag) then
3451              i01=ii
3452              liiflag=.true.
3453           endif
3454          enddo
3455         enddo
3456         lim_odl=iistart-1
3457       endif
3458 c      write (iout,*) "Removing distances completed"
3459 c      call flush(iout)
3460       endif ! .not. klapaucjusz
3461
3462       if (constr_homology.gt.0) call homology_partition
3463 c      write (iout,*) "After homology_partition"
3464 c      call flush(iout)
3465       if (constr_homology.gt.0) call init_int_table
3466 c      write (iout,*) "After init_int_table"
3467 c      call flush(iout)
3468 c      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3469 c      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3470 c
3471 c Print restraints
3472 c
3473       if (.not.out_template_restr) return
3474 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3475       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3476        write (iout,*) "Distance restraints from templates"
3477        do ii=1,lim_odl
3478        write(iout,'(3i7,100(2f8.2,1x,l1,4x))')
3479      &  ii,ires_homo(ii),jres_homo(ii),
3480      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3481      &  ki=1,constr_homology)
3482        enddo
3483        write (iout,*) "Dihedral angle restraints from templates"
3484        do i=nnt+3,nct
3485         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3486      &      (rad2deg*dih(ki,i),
3487      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3488        enddo
3489        write (iout,*) "Virtual-bond angle restraints from templates"
3490        do i=nnt+2,nct
3491         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3492      &      (rad2deg*thetatpl(ki,i),
3493      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3494        enddo
3495        write (iout,*) "SC restraints from templates"
3496        do i=nnt,nct
3497         write(iout,'(i7,100(4f8.2,4x))') i,
3498      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3499      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3500        enddo
3501       endif
3502 c -----------------------------------------------------------------
3503       return
3504       end
3505 c----------------------------------------------------------------------
3506 #ifdef WINIFL
3507       subroutine flush(iu)
3508       return
3509       end
3510 #endif
3511 #ifdef AIX
3512       subroutine flush(iu)
3513       call flush_(iu)
3514       return
3515       end
3516 #endif
3517 c------------------------------------------------------------------------------
3518       subroutine copy_to_tmp(source)
3519       implicit none
3520       include "DIMENSIONS"
3521       include "COMMON.IOUNITS"
3522       character*(*) source
3523       character* 256 tmpfile
3524       integer ilen
3525       external ilen
3526       logical ex
3527       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3528       inquire(file=tmpfile,exist=ex)
3529       if (ex) then
3530         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3531      &   " to temporary directory..."
3532         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3533         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3534       endif
3535       return
3536       end
3537 c------------------------------------------------------------------------------
3538       subroutine move_from_tmp(source)
3539       implicit none
3540       include "DIMENSIONS"
3541       include "COMMON.IOUNITS"
3542       character*(*) source
3543       integer ilen
3544       external ilen
3545       write (*,*) "Moving ",source(:ilen(source)),
3546      & " from temporary directory to working directory"
3547       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3548       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3549       return
3550       end
3551 c------------------------------------------------------------------------------
3552       subroutine random_init(seed)
3553 C
3554 C Initialize random number generator
3555 C
3556       implicit none
3557       include 'DIMENSIONS'
3558 #ifdef MPI
3559       include 'mpif.h'
3560       logical OKRandom, prng_restart
3561       real*8  r1
3562       integer iseed_array(4)
3563       integer error_msg,ierr
3564 #endif
3565       include 'COMMON.IOUNITS'
3566       include 'COMMON.TIME1'
3567       include 'COMMON.THREAD'
3568       include 'COMMON.SBRIDGE'
3569       include 'COMMON.CONTROL'
3570       include 'COMMON.MCM'
3571       include 'COMMON.MAP'
3572       include 'COMMON.HEADER'
3573       include 'COMMON.CSA'
3574       include 'COMMON.CHAIN'
3575       include 'COMMON.MUCA'
3576       include 'COMMON.MD'
3577       include 'COMMON.FFIELD'
3578       include 'COMMON.SETUP'
3579       integer i,iseed
3580       double precision seed,ran_number
3581       iseed=-dint(dabs(seed))
3582       if (iseed.eq.0) then
3583         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3584      &    'Random seed undefined. The program will stop.'
3585         write (*,'(/80(1h*)/20x,a/80(1h*))')
3586      &    'Random seed undefined. The program will stop.'
3587 #ifdef MPI
3588         call mpi_finalize(mpi_comm_world,ierr)
3589 #endif
3590         stop 'Bad random seed.'
3591       endif
3592 #ifdef MPI
3593       if (fg_rank.eq.0) then
3594       seed=seed*(me+1)+1
3595 #ifdef AMD64
3596       if(me.eq.king)
3597      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3598       OKRandom = prng_restart(me,iseed)
3599 #else
3600       do i=1,4
3601        tmp=65536.0d0**(4-i)
3602        iseed_array(i) = dint(seed/tmp)
3603        seed=seed-iseed_array(i)*tmp
3604       enddo
3605       if(me.eq.king)
3606      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3607      &                 (iseed_array(i),i=1,4)
3608       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3609      &                 (iseed_array(i),i=1,4)
3610       OKRandom = prng_restart(me,iseed_array)
3611 #endif
3612       if (OKRandom) then
3613 c        r1 = prng_next(me)
3614         r1=ran_number(0.0D0,1.0D0)
3615         if(me.eq.king)
3616      &   write (iout,*) 'ran_num',r1
3617         if (r1.lt.0.0d0) OKRandom=.false.
3618       endif
3619       if (.not.OKRandom) then
3620         write (iout,*) 'PRNG IS NOT WORKING!!!'
3621         print *,'PRNG IS NOT WORKING!!!'
3622         if (me.eq.0) then
3623          call flush(iout)
3624          call mpi_abort(mpi_comm_world,error_msg,ierr)
3625          stop
3626         else
3627          write (iout,*) 'too many processors for parallel prng'
3628          write (*,*) 'too many processors for parallel prng'
3629          call flush(iout)
3630          stop
3631         endif
3632       endif
3633       endif
3634 #else
3635       call vrndst(iseed)
3636       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3637 #endif
3638       return
3639       end
3640 c----------------------------------------------------------------------
3641       subroutine read_klapaucjusz
3642       implicit none
3643       include 'DIMENSIONS'
3644 #ifdef MPI
3645       include 'mpif.h'
3646 #endif
3647       include 'COMMON.SETUP'
3648       include 'COMMON.CONTROL'
3649       include 'COMMON.HOMOLOGY'
3650       include 'COMMON.CHAIN'
3651       include 'COMMON.IOUNITS'
3652       include 'COMMON.MD'
3653       include 'COMMON.GEO'
3654       include 'COMMON.INTERACT'
3655       include 'COMMON.NAMES'
3656       character*256 fragfile
3657       integer ninclust(maxclust),inclust(max_template,maxclust),
3658      &  nresclust(maxclust),iresclust(maxres,maxclust),nclust
3659
3660       character*2 kic2
3661       character*24 model_ki_dist, model_ki_angle
3662       character*500 controlcard
3663       integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp,
3664      & idomain_tmp,
3665      & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,
3666      & i01,i10,nnt_chain,nct_chain
3667       integer itype_temp(maxres)
3668       double precision distal
3669       logical lprn /.true./
3670       integer nres_temp
3671       integer ilen
3672       external ilen
3673       logical liiflag,lfirst
3674 c
3675 c
3676       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3677       double precision, dimension (max_template,maxres) :: rescore
3678       double precision, dimension (max_template,maxres) :: rescore2
3679       character*24 tpl_k_rescore
3680       character*256 pdbfile
3681
3682 c
3683 c For new homol impl
3684 c
3685       include 'COMMON.VAR'
3686 c
3687 c      write (iout,*) "READ_KLAPAUCJUSZ"
3688 c      print *,"READ_KLAPAUCJUSZ"
3689 c      call flush(iout)
3690       call getenv("FRAGFILE",fragfile)
3691       write (iout,*) "Opening", fragfile
3692       call flush(iout)
3693       open(ientin,file=fragfile,status="old",err=10)
3694 c      write (iout,*) " opened"
3695 c      call flush(iout)
3696
3697       sigma_theta=0.0
3698       sigma_d=0.0
3699       sigma_dih=0.0
3700       l_homo = .false.
3701
3702       nres_temp=nres
3703       itype_temp=itype
3704       ii=0
3705       lim_odl=0
3706
3707 c      write (iout,*) "Entering loop"
3708 c      call flush(iout)
3709
3710       DO IGR = 1,NCHAIN_GROUP
3711
3712 c      write (iout,*) "igr",igr
3713       call flush(iout)
3714       read(ientin,*) constr_homology,nclust
3715
3716       if (start_from_model) then
3717         nmodel_start=constr_homology
3718       else
3719         nmodel_start=0
3720       endif
3721
3722       ii_old=lim_odl
3723
3724       ichain=iequiv(1,igr)
3725       nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
3726       nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
3727 c      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
3728
3729 c Read pdb files
3730       do k=1,constr_homology
3731         read(ientin,'(a)') pdbfile
3732         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3733      &  pdbfile(:ilen(pdbfile))
3734         pdbfiles_chomo(k)=pdbfile
3735         open(ipdbin,file=pdbfile,status='old',err=33)
3736         goto 34
3737   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3738      &  pdbfile(:ilen(pdbfile))
3739         stop
3740   34    continue
3741         unres_pdb=.false.
3742 c        nres_temp=nres
3743         call readpdb_template(k)
3744         nres_chomo(k)=nres
3745 c        nres=nres_temp
3746         do i=1,nres
3747           rescore(k,i)=0.2d0
3748           rescore2(k,i)=1.0d0
3749         enddo
3750       enddo
3751 c Read clusters
3752       do i=1,nclust
3753         read(ientin,*) ninclust(i),nresclust(i)
3754         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3755         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3756       enddo
3757 c
3758 c Loop over clusters
3759 c
3760       do l=1,nclust
3761         do ll = 1,ninclust(l)
3762
3763         k = inclust(ll,l)
3764 c        write (iout,*) "l",l," ll",ll," k",k
3765         do i=1,nres
3766           idomain(k,i)=0
3767         enddo
3768         do i=1,nresclust(l)
3769           if (nnt.gt.1)  then
3770             idomain(k,iresclust(i,l)+1) = 1
3771           else
3772             idomain(k,iresclust(i,l)) = 1
3773           endif
3774         enddo
3775 c
3776 c     Distance restraints
3777 c
3778 c          ... --> odl(k,ii)
3779 C Copy the coordinates from reference coordinates (?)
3780 c        nres_temp=nres
3781         nres=nres_chomo(k)
3782         do i=1,2*nres
3783           do j=1,3
3784             c(j,i)=chomo(j,i,k)
3785 c           write (iout,*) "c(",j,i,") =",c(j,i)
3786           enddo
3787         enddo
3788         call int_from_cart(.true.,.false.)
3789         call sc_loc_geom(.false.)
3790         do i=1,nres
3791           thetaref(i)=theta(i)
3792           phiref(i)=phi(i)
3793         enddo
3794 c        nres=nres_temp
3795         if (waga_dist.ne.0.0d0) then
3796           ii=ii_old
3797 c          do i = nnt,nct-2
3798           do i = nnt_chain,nct_chain-2
3799 c            do j=i+2,nct
3800             do j=i+2,nct_chain
3801
3802               x12=c(1,i)-c(1,j)
3803               y12=c(2,i)-c(2,j)
3804               z12=c(3,i)-c(3,j)
3805               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3806 c              write (iout,*) k,i,j,distal,dist2_cut
3807
3808             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3809      &            .and. distal.le.dist2_cut ) then
3810
3811               ii=ii+1
3812               ii_in_use(ii)=1
3813               l_homo(k,ii)=.true.
3814
3815 c             write (iout,*) "k",k
3816 c             write (iout,*) "i",i," j",j," constr_homology",
3817 c     &                       constr_homology
3818               ires_homo(ii)=i+chain_border1(1,igr)-1
3819               jres_homo(ii)=j+chain_border1(1,igr)-1
3820               odl(k,ii)=distal
3821               if (read2sigma) then
3822                 sigma_odl(k,ii)=0
3823                 do ik=i,j
3824                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3825                 enddo
3826                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3827                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3828      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3829               else
3830                 if (odl(k,ii).le.dist_cut) then
3831                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3832                 else
3833 #ifdef OLDSIGMA
3834                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3835      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3836 #else
3837                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3838      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3839 #endif
3840                 endif
3841               endif
3842               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3843             else
3844               ii=ii+1
3845 c              l_homo(k,ii)=.false.
3846             endif
3847             enddo
3848           enddo
3849         lim_odl=ii
3850         endif
3851 c
3852 c     Theta, dihedral and SC retraints
3853 c
3854         if (waga_angle.gt.0.0d0) then
3855           do i = nnt_chain+3,nct_chain
3856             iii=i+chain_border1(1,igr)-1
3857             if (idomain(k,i).eq.0) then
3858 c               sigma_dih(k,i)=0.0
3859                cycle
3860             endif
3861             dih(k,iii)=phiref(i)
3862             sigma_dih(k,iii)=
3863      &         (rescore(k,i)+rescore(k,i-1)+
3864      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3865 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3866 c     &       " sigma_dihed",sigma_dih(k,i)
3867             if (sigma_dih(k,iii).ne.0)
3868      &       sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
3869           enddo
3870 c          lim_dih=nct-nnt-2
3871         endif
3872
3873         if (waga_theta.gt.0.0d0) then
3874           do i = nnt_chain+2,nct_chain
3875              iii=i+chain_border1(1,igr)-1
3876              if (idomain(k,i).eq.0) then
3877 c              sigma_theta(k,i)=0.0
3878               cycle
3879              endif
3880              thetatpl(k,iii)=thetaref(i)
3881              sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+
3882      &                        rescore(k,i-2))/3.0
3883              if (sigma_theta(k,iii).ne.0)
3884      &       sigma_theta(k,iii)=1.0d0/
3885      &       (sigma_theta(k,iii)*sigma_theta(k,iii))
3886           enddo
3887         endif
3888
3889         if (waga_d.gt.0.0d0) then
3890           do i = nnt_chain,nct_chain
3891              iii=i+chain_border1(1,igr)-1
3892                if (itype(i).eq.10) cycle
3893                if (idomain(k,i).eq.0 ) then
3894 c                  sigma_d(k,i)=0.0
3895                   cycle
3896                endif
3897                xxtpl(k,iii)=xxref(i)
3898                yytpl(k,iii)=yyref(i)
3899                zztpl(k,iii)=zzref(i)
3900                sigma_d(k,iii)=rescore(k,i)
3901                if (sigma_d(k,iii).ne.0)
3902      &          sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
3903 c               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3904           enddo
3905         endif
3906       enddo ! l
3907       enddo ! ll
3908 c
3909 c remove distance restraints not used in any model from the list
3910 c shift data in all arrays
3911 c
3912 c      write (iout,*) "ii_old",ii_old
3913       if (waga_dist.ne.0.0d0) then
3914 #ifdef DEBUG
3915        write (iout,*) "Distance restraints from templates"
3916        do iii=1,lim_odl
3917        write(iout,'(4i5,100(2f8.2,1x,l1,4x))')
3918      &  iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii),
3919      &  (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii),
3920      &  ki=1,constr_homology)
3921        enddo
3922 #endif
3923         ii=ii_old
3924         liiflag=.true.
3925         lfirst=.true.
3926         do i=nnt_chain,nct_chain-2
3927          do j=i+2,nct_chain
3928           ii=ii+1
3929 c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3930 c     &            .and. distal.le.dist2_cut ) then
3931 c          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
3932 c          call flush(iout)
3933           if (ii_in_use(ii).eq.0.and.liiflag.or.
3934      &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3935             liiflag=.false.
3936             i10=ii
3937             if (lfirst) then
3938               lfirst=.false.
3939               iistart=ii
3940             else
3941               if(i10.eq.lim_odl) i10=i10+1
3942               do ki=0,i10-i01-1
3943                ires_homo(iistart+ki)=ires_homo(ki+i01)
3944                jres_homo(iistart+ki)=jres_homo(ki+i01)
3945                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3946                do k=1,constr_homology
3947                 odl(k,iistart+ki)=odl(k,ki+i01)
3948                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3949                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3950                enddo
3951               enddo
3952               iistart=iistart+i10-i01
3953             endif
3954           endif
3955           if (ii_in_use(ii).ne.0.and..not.liiflag) then
3956              i01=ii
3957              liiflag=.true.
3958           endif
3959          enddo
3960         enddo
3961         lim_odl=iistart-1
3962       endif
3963
3964       lll=lim_odl-ii_old
3965
3966       do i=2,nequiv(igr)
3967
3968         ichain=iequiv(i,igr)
3969
3970         do j=nnt_chain,nct_chain
3971           jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
3972           do k=1,constr_homology
3973             dih(k,jj)=dih(k,j)
3974             sigma_dih(k,jj)=sigma_dih(k,j)
3975             thetatpl(k,jj)=thetatpl(k,j)
3976             sigma_theta(k,jj)=sigma_theta(k,j)
3977             xxtpl(k,jj)=xxtpl(k,j)
3978             yytpl(k,jj)=yytpl(k,j)
3979             zztpl(k,jj)=zztpl(k,j)
3980             sigma_d(k,jj)=sigma_d(k,j)
3981           enddo
3982         enddo
3983
3984         jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
3985 c        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
3986         do j=ii_old+1,lim_odl
3987           ires_homo(j+lll)=ires_homo(j)+jj
3988           jres_homo(j+lll)=jres_homo(j)+jj
3989           do k=1,constr_homology
3990             odl(k,j+lll)=odl(k,j)
3991             sigma_odl(k,j+lll)=sigma_odl(k,j)
3992             l_homo(k,j+lll)=l_homo(k,j)
3993           enddo
3994         enddo
3995
3996         ii_old=ii_old+lll
3997         lim_odl=lim_odl+lll
3998
3999       enddo
4000
4001       ENDDO ! IGR
4002
4003       if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
4004       nres=nres_temp
4005       itype=itype_temp
4006
4007       return
4008    10 stop "Error in fragment file"
4009       end