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       dist1cut=(index(controlcard,'DIST1CUT').gt.0)
3109       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3110       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3111       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3112       if(.not.read2sigma.and.start_from_model) then
3113           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3114      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3115           start_from_model=.false.
3116           iranconf=(indpdb.le.0)
3117       endif
3118       if(start_from_model .and. (me.eq.king .or. .not. out1file))
3119      &    write(iout,*) 'START_FROM_MODELS is ON'
3120 c      if(start_from_model .and. rest) then
3121 c        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3122 c         write(iout,*) 'START_FROM_MODELS is OFF'
3123 c         write(iout,*) 'remove restart keyword from input'
3124 c        endif
3125 c      endif
3126       if (start_from_model) nmodel_start=constr_homology
3127       if (homol_nset.gt.1)then
3128          call card_concat(controlcard)
3129          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3130          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3131 c          write(iout,*) "iset homology_weight "
3132           do i=1,homol_nset
3133            write(iout,*) i,waga_homology(i)
3134           enddo
3135          endif
3136          iset=mod(kolor,homol_nset)+1
3137       else
3138        iset=1
3139        waga_homology(1)=1.0
3140       endif
3141
3142 cd      write (iout,*) "nnt",nnt," nct",nct
3143 cd      call flush(iout)
3144
3145       if (read_homol_frag) then
3146         call read_klapaucjusz
3147       else
3148
3149       lim_odl=0
3150       lim_dih=0
3151 c
3152 c      write(iout,*) 'nnt=',nnt,'nct=',nct
3153 c
3154 c      do i = nnt,nct
3155 c        do k=1,constr_homology
3156 c         idomain(k,i)=0
3157 c        enddo
3158 c      enddo
3159        idomain=0
3160
3161 c      ii=0
3162 c      do i = nnt,nct-2
3163 c        do j=i+2,nct
3164 c        ii=ii+1
3165 c        ii_in_use(ii)=0
3166 c        enddo
3167 c      enddo
3168       ii_in_use=0
3169       do k=1,constr_homology
3170
3171         read(inp,'(a)') pdbfile
3172         pdbfiles_chomo(k)=pdbfile
3173         if(me.eq.king .or. .not. out1file)
3174      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3175      &  pdbfile(:ilen(pdbfile))
3176         open(ipdbin,file=pdbfile,status='old',err=33)
3177         goto 34
3178   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3179      &  pdbfile(:ilen(pdbfile))
3180         stop
3181   34    continue
3182 c        print *,'Begin reading pdb data'
3183 c
3184 c Files containing res sim or local scores (former containing sigmas)
3185 c
3186
3187         write(kic2,'(bz,i2.2)') k
3188
3189         tpl_k_rescore="template"//kic2//".sco"
3190
3191         unres_pdb=.false.
3192         nres_temp=nres
3193         if (read2sigma) then
3194           call readpdb_template(k)
3195         else
3196           call readpdb
3197         endif
3198         nres_chomo(k)=nres
3199         nres=nres_temp
3200 c
3201 c     Distance restraints
3202 c
3203 c          ... --> odl(k,ii)
3204 C Copy the coordinates from reference coordinates (?)
3205         do i=1,2*nres_chomo(k)
3206           do j=1,3
3207             c(j,i)=cref(j,i)
3208 c           write (iout,*) "c(",j,i,") =",c(j,i)
3209           enddo
3210         enddo
3211 c
3212 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3213 c
3214 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3215           open (ientin,file=tpl_k_rescore,status='old')
3216           if (nnt.gt.1) rescore(k,1)=0.0d0
3217           do irec=nnt,nct ! loop for reading res sim
3218             if (read2sigma) then
3219              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3220      &                                rescore3_tmp,idomain_tmp
3221              i_tmp=i_tmp+nnt-1
3222              idomain(k,i_tmp)=idomain_tmp
3223              rescore(k,i_tmp)=rescore_tmp
3224              rescore2(k,i_tmp)=rescore2_tmp
3225              rescore3(k,i_tmp)=rescore3_tmp
3226              if (.not. out1file .or. me.eq.king)
3227      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3228      &                      i_tmp,rescore2_tmp,rescore_tmp,
3229      &                                rescore3_tmp,idomain_tmp
3230             else
3231              idomain(k,irec)=1
3232              read (ientin,*,end=1401) rescore_tmp
3233
3234 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3235              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3236 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3237             endif
3238           enddo
3239  1401   continue
3240         close (ientin)
3241         if (waga_dist.ne.0.0d0) then
3242           ii=0
3243           do i = nnt,nct-2
3244             do j=i+2,nct
3245
3246               x12=c(1,i)-c(1,j)
3247               y12=c(2,i)-c(2,j)
3248               z12=c(3,i)-c(3,j)
3249               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3250 c              write (iout,*) k,i,j,distal,dist2_cut
3251            if (dist1cut .and. k.gt.1) then
3252               ii=ii+1
3253               if (l_homo(1,ii)) then
3254                 ii_in_use(ii)=1
3255                 l_homo(k,ii)=.true.
3256                 ires_homo(ii)=i
3257                 jres_homo(ii)=j
3258                 odl(k,ii)=distal
3259                 sigma_odl(k,ii)=sigma_odl(1,ii)
3260               else
3261                 l_homo(k,ii)=.false.
3262               endif
3263            else
3264             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3265      &            .and. distal.le.dist2_cut ) then
3266
3267               ii=ii+1
3268               ii_in_use(ii)=1
3269               l_homo(k,ii)=.true.
3270
3271 c             write (iout,*) "k",k
3272 c             write (iout,*) "i",i," j",j," constr_homology",
3273 c    &                       constr_homology
3274               ires_homo(ii)=i
3275               jres_homo(ii)=j
3276               odl(k,ii)=distal
3277               if (read2sigma) then
3278                 sigma_odl(k,ii)=0
3279                 do ik=i,j
3280                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3281                 enddo
3282                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3283                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3284      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3285               else
3286                 if (odl(k,ii).le.dist_cut) then
3287                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3288                 else
3289 #ifdef OLDSIGMA
3290                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3291      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3292 #else
3293                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3294      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3295 #endif
3296                 endif
3297               endif
3298               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3299             else
3300 c              ii=ii+1
3301 c              l_homo(k,ii)=.false.
3302             endif
3303            endif
3304             enddo
3305           enddo
3306         lim_odl=ii
3307         endif
3308 c        write (iout,*) "Distance restraints set"
3309 c        call flush(iout)
3310 c
3311 c     Theta, dihedral and SC retraints
3312 c
3313         if (waga_angle.gt.0.0d0) then
3314 c         open (ientin,file=tpl_k_sigma_dih,status='old')
3315 c         do irec=1,maxres-3 ! loop for reading sigma_dih
3316 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3317 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3318 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3319 c    &                            sigma_dih(k,i+nnt-1)
3320 c         enddo
3321 c1402   continue
3322 c         close (ientin)
3323           do i = nnt+3,nct
3324             if (idomain(k,i).eq.0) then
3325                sigma_dih(k,i)=0.0
3326                cycle
3327             endif
3328             dih(k,i)=phiref(i) ! right?
3329 c           read (ientin,*) sigma_dih(k,i) ! original variant
3330 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
3331 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3332 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3333 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
3334 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
3335
3336             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3337      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3338 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3339 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
3340 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3341 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3342 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
3343 c   Instead of res sim other local measure of b/b str reliability possible
3344             if (sigma_dih(k,i).ne.0)
3345      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3346 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3347           enddo
3348           lim_dih=nct-nnt-2
3349         endif
3350 c        write (iout,*) "Dihedral angle restraints set"
3351 c        call flush(iout)
3352
3353         if (waga_theta.gt.0.0d0) then
3354 c         open (ientin,file=tpl_k_sigma_theta,status='old')
3355 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3356 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3357 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3358 c    &                              sigma_theta(k,i+nnt-1)
3359 c         enddo
3360 c1403   continue
3361 c         close (ientin)
3362
3363           do i = nnt+2,nct ! right? without parallel.
3364 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
3365 c         do i=ithet_start,ithet_end ! with FG parallel.
3366              if (idomain(k,i).eq.0) then
3367               sigma_theta(k,i)=0.0
3368               cycle
3369              endif
3370              thetatpl(k,i)=thetaref(i)
3371 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3372 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
3373 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3374 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
3375 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
3376              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3377      &                        rescore(k,i-2))/3.0
3378 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3379              if (sigma_theta(k,i).ne.0)
3380      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3381
3382 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3383 c                             rescore(k,i-2) !  right expression ?
3384 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3385           enddo
3386         endif
3387 c        write (iout,*) "Angle restraints set"
3388 c        call flush(iout)
3389
3390         if (waga_d.gt.0.0d0) then
3391 c       open (ientin,file=tpl_k_sigma_d,status='old')
3392 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3393 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3394 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3395 c    &                          sigma_d(k,i+nnt-1)
3396 c         enddo
3397 c1404   continue
3398
3399           do i = nnt,nct ! right? without parallel.
3400 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3401 c         do i=loc_start,loc_end ! with FG parallel.
3402                if (itype(i).eq.10) cycle
3403                if (idomain(k,i).eq.0 ) then
3404                   sigma_d(k,i)=0.0
3405                   cycle
3406                endif
3407                xxtpl(k,i)=xxref(i)
3408                yytpl(k,i)=yyref(i)
3409                zztpl(k,i)=zzref(i)
3410 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3411 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3412 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3413 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3414                sigma_d(k,i)=rescore3(k,i) !  right expression ?
3415                if (sigma_d(k,i).ne.0)
3416      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3417
3418 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3419 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3420 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3421           enddo
3422         endif
3423       enddo
3424 c      write (iout,*) "SC restraints set"
3425 c      call flush(iout)
3426 c
3427 c remove distance restraints not used in any model from the list
3428 c shift data in all arrays
3429 c
3430 c      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3431       if (waga_dist.ne.0.0d0) then
3432         ii=0
3433         liiflag=.true.
3434         lfirst=.true.
3435         do i=nnt,nct-2
3436          do j=i+2,nct
3437           ii=ii+1
3438 c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3439 c     &            .and. distal.le.dist2_cut ) then
3440 c          write (iout,*) "i",i," j",j," ii",ii
3441 c          call flush(iout)
3442           if (ii_in_use(ii).eq.0.and.liiflag.or.
3443      &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3444             liiflag=.false.
3445             i10=ii
3446             if (lfirst) then
3447               lfirst=.false.
3448               iistart=ii
3449             else
3450               if(i10.eq.lim_odl) i10=i10+1
3451               do ki=0,i10-i01-1
3452                ires_homo(iistart+ki)=ires_homo(ki+i01)
3453                jres_homo(iistart+ki)=jres_homo(ki+i01)
3454                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3455                do k=1,constr_homology
3456                 odl(k,iistart+ki)=odl(k,ki+i01)
3457                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3458                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3459                enddo
3460               enddo
3461               iistart=iistart+i10-i01
3462             endif
3463           endif
3464           if (ii_in_use(ii).ne.0.and..not.liiflag) then
3465              i01=ii
3466              liiflag=.true.
3467           endif
3468          enddo
3469         enddo
3470         lim_odl=iistart-1
3471       endif
3472 c      write (iout,*) "Removing distances completed"
3473 c      call flush(iout)
3474       endif ! .not. klapaucjusz
3475
3476       if (constr_homology.gt.0) call homology_partition
3477 c      write (iout,*) "After homology_partition"
3478 c      call flush(iout)
3479       if (constr_homology.gt.0) call init_int_table
3480 c      write (iout,*) "After init_int_table"
3481 c      call flush(iout)
3482 c      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3483 c      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3484 c
3485 c Print restraints
3486 c
3487       if (.not.out_template_restr) return
3488 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3489       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3490        write (iout,*) "Distance restraints from templates"
3491        do ii=1,lim_odl
3492        write(iout,'(3i7,100(2f8.2,1x,l1,4x))')
3493      &  ii,ires_homo(ii),jres_homo(ii),
3494      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3495      &  ki=1,constr_homology)
3496        enddo
3497        write (iout,*) "Dihedral angle restraints from templates"
3498        do i=nnt+3,nct
3499         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3500      &      (rad2deg*dih(ki,i),
3501      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3502        enddo
3503        write (iout,*) "Virtual-bond angle restraints from templates"
3504        do i=nnt+2,nct
3505         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3506      &      (rad2deg*thetatpl(ki,i),
3507      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3508        enddo
3509        write (iout,*) "SC restraints from templates"
3510        do i=nnt,nct
3511         write(iout,'(i7,100(4f8.2,4x))') i,
3512      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3513      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3514        enddo
3515       endif
3516 c -----------------------------------------------------------------
3517       return
3518       end
3519 c----------------------------------------------------------------------
3520 #ifdef WINIFL
3521       subroutine flush(iu)
3522       return
3523       end
3524 #endif
3525 #ifdef AIX
3526       subroutine flush(iu)
3527       call flush_(iu)
3528       return
3529       end
3530 #endif
3531 c------------------------------------------------------------------------------
3532       subroutine copy_to_tmp(source)
3533       implicit none
3534       include "DIMENSIONS"
3535       include "COMMON.IOUNITS"
3536       character*(*) source
3537       character* 256 tmpfile
3538       integer ilen
3539       external ilen
3540       logical ex
3541       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3542       inquire(file=tmpfile,exist=ex)
3543       if (ex) then
3544         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3545      &   " to temporary directory..."
3546         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3547         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3548       endif
3549       return
3550       end
3551 c------------------------------------------------------------------------------
3552       subroutine move_from_tmp(source)
3553       implicit none
3554       include "DIMENSIONS"
3555       include "COMMON.IOUNITS"
3556       character*(*) source
3557       integer ilen
3558       external ilen
3559       write (*,*) "Moving ",source(:ilen(source)),
3560      & " from temporary directory to working directory"
3561       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3562       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3563       return
3564       end
3565 c------------------------------------------------------------------------------
3566       subroutine random_init(seed)
3567 C
3568 C Initialize random number generator
3569 C
3570       implicit none
3571       include 'DIMENSIONS'
3572 #ifdef MPI
3573       include 'mpif.h'
3574       logical OKRandom, prng_restart
3575       real*8  r1
3576       integer iseed_array(4)
3577       integer error_msg,ierr
3578 #endif
3579       include 'COMMON.IOUNITS'
3580       include 'COMMON.TIME1'
3581       include 'COMMON.THREAD'
3582       include 'COMMON.SBRIDGE'
3583       include 'COMMON.CONTROL'
3584       include 'COMMON.MCM'
3585       include 'COMMON.MAP'
3586       include 'COMMON.HEADER'
3587       include 'COMMON.CSA'
3588       include 'COMMON.CHAIN'
3589       include 'COMMON.MUCA'
3590       include 'COMMON.MD'
3591       include 'COMMON.FFIELD'
3592       include 'COMMON.SETUP'
3593       integer i,iseed
3594       double precision seed,ran_number
3595       iseed=-dint(dabs(seed))
3596       if (iseed.eq.0) then
3597         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3598      &    'Random seed undefined. The program will stop.'
3599         write (*,'(/80(1h*)/20x,a/80(1h*))')
3600      &    'Random seed undefined. The program will stop.'
3601 #ifdef MPI
3602         call mpi_finalize(mpi_comm_world,ierr)
3603 #endif
3604         stop 'Bad random seed.'
3605       endif
3606 #ifdef MPI
3607       if (fg_rank.eq.0) then
3608       seed=seed*(me+1)+1
3609 #ifdef AMD64
3610       if(me.eq.king)
3611      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3612       OKRandom = prng_restart(me,iseed)
3613 #else
3614       do i=1,4
3615        tmp=65536.0d0**(4-i)
3616        iseed_array(i) = dint(seed/tmp)
3617        seed=seed-iseed_array(i)*tmp
3618       enddo
3619       if(me.eq.king)
3620      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3621      &                 (iseed_array(i),i=1,4)
3622       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3623      &                 (iseed_array(i),i=1,4)
3624       OKRandom = prng_restart(me,iseed_array)
3625 #endif
3626       if (OKRandom) then
3627 c        r1 = prng_next(me)
3628         r1=ran_number(0.0D0,1.0D0)
3629         if(me.eq.king)
3630      &   write (iout,*) 'ran_num',r1
3631         if (r1.lt.0.0d0) OKRandom=.false.
3632       endif
3633       if (.not.OKRandom) then
3634         write (iout,*) 'PRNG IS NOT WORKING!!!'
3635         print *,'PRNG IS NOT WORKING!!!'
3636         if (me.eq.0) then
3637          call flush(iout)
3638          call mpi_abort(mpi_comm_world,error_msg,ierr)
3639          stop
3640         else
3641          write (iout,*) 'too many processors for parallel prng'
3642          write (*,*) 'too many processors for parallel prng'
3643          call flush(iout)
3644          stop
3645         endif
3646       endif
3647       endif
3648 #else
3649       call vrndst(iseed)
3650       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3651 #endif
3652       return
3653       end
3654 c----------------------------------------------------------------------
3655       subroutine read_klapaucjusz
3656       implicit none
3657       include 'DIMENSIONS'
3658 #ifdef MPI
3659       include 'mpif.h'
3660 #endif
3661       include 'COMMON.SETUP'
3662       include 'COMMON.CONTROL'
3663       include 'COMMON.HOMOLOGY'
3664       include 'COMMON.CHAIN'
3665       include 'COMMON.IOUNITS'
3666       include 'COMMON.MD'
3667       include 'COMMON.GEO'
3668       include 'COMMON.INTERACT'
3669       include 'COMMON.NAMES'
3670       character*256 fragfile
3671       integer ninclust(maxclust),inclust(max_template,maxclust),
3672      &  nresclust(maxclust),iresclust(maxres,maxclust),nclust
3673
3674       character*2 kic2
3675       character*24 model_ki_dist, model_ki_angle
3676       character*500 controlcard
3677       integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp,
3678      & idomain_tmp,
3679      & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,
3680      & i01,i10,nnt_chain,nct_chain
3681       integer itype_temp(maxres)
3682       double precision distal
3683       logical lprn /.true./
3684       integer nres_temp
3685       integer ilen
3686       external ilen
3687       logical liiflag,lfirst
3688 c
3689 c
3690       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3691       double precision, dimension (max_template,maxres) :: rescore
3692       double precision, dimension (max_template,maxres) :: rescore2
3693       character*24 tpl_k_rescore
3694       character*256 pdbfile
3695
3696 c
3697 c For new homol impl
3698 c
3699       include 'COMMON.VAR'
3700 c
3701 c      write (iout,*) "READ_KLAPAUCJUSZ"
3702 c      print *,"READ_KLAPAUCJUSZ"
3703 c      call flush(iout)
3704       call getenv("FRAGFILE",fragfile)
3705       write (iout,*) "Opening", fragfile
3706       call flush(iout)
3707       open(ientin,file=fragfile,status="old",err=10)
3708 c      write (iout,*) " opened"
3709 c      call flush(iout)
3710
3711       sigma_theta=0.0
3712       sigma_d=0.0
3713       sigma_dih=0.0
3714       l_homo = .false.
3715
3716       nres_temp=nres
3717       itype_temp=itype
3718       ii=0
3719       lim_odl=0
3720
3721 c      write (iout,*) "Entering loop"
3722 c      call flush(iout)
3723
3724       DO IGR = 1,NCHAIN_GROUP
3725
3726 c      write (iout,*) "igr",igr
3727       call flush(iout)
3728       read(ientin,*) constr_homology,nclust
3729
3730       if (start_from_model) then
3731         nmodel_start=constr_homology
3732       else
3733         nmodel_start=0
3734       endif
3735
3736       ii_old=lim_odl
3737
3738       ichain=iequiv(1,igr)
3739       nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
3740       nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
3741 c      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
3742
3743 c Read pdb files
3744       do k=1,constr_homology
3745         read(ientin,'(a)') pdbfile
3746         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3747      &  pdbfile(:ilen(pdbfile))
3748         pdbfiles_chomo(k)=pdbfile
3749         open(ipdbin,file=pdbfile,status='old',err=33)
3750         goto 34
3751   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3752      &  pdbfile(:ilen(pdbfile))
3753         stop
3754   34    continue
3755         unres_pdb=.false.
3756 c        nres_temp=nres
3757         call readpdb_template(k)
3758         nres_chomo(k)=nres
3759 c        nres=nres_temp
3760         do i=1,nres
3761           rescore(k,i)=0.2d0
3762           rescore2(k,i)=1.0d0
3763         enddo
3764       enddo
3765 c Read clusters
3766       do i=1,nclust
3767         read(ientin,*) ninclust(i),nresclust(i)
3768         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3769         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3770       enddo
3771 c
3772 c Loop over clusters
3773 c
3774       do l=1,nclust
3775         do ll = 1,ninclust(l)
3776
3777         k = inclust(ll,l)
3778 c        write (iout,*) "l",l," ll",ll," k",k
3779         do i=1,nres
3780           idomain(k,i)=0
3781         enddo
3782         do i=1,nresclust(l)
3783           if (nnt.gt.1)  then
3784             idomain(k,iresclust(i,l)+1) = 1
3785           else
3786             idomain(k,iresclust(i,l)) = 1
3787           endif
3788         enddo
3789 c
3790 c     Distance restraints
3791 c
3792 c          ... --> odl(k,ii)
3793 C Copy the coordinates from reference coordinates (?)
3794 c        nres_temp=nres
3795         nres=nres_chomo(k)
3796         do i=1,2*nres
3797           do j=1,3
3798             c(j,i)=chomo(j,i,k)
3799 c           write (iout,*) "c(",j,i,") =",c(j,i)
3800           enddo
3801         enddo
3802         call int_from_cart(.true.,.false.)
3803         call sc_loc_geom(.false.)
3804         do i=1,nres
3805           thetaref(i)=theta(i)
3806           phiref(i)=phi(i)
3807         enddo
3808 c        nres=nres_temp
3809         if (waga_dist.ne.0.0d0) then
3810           ii=ii_old
3811 c          do i = nnt,nct-2
3812           do i = nnt_chain,nct_chain-2
3813 c            do j=i+2,nct
3814             do j=i+2,nct_chain
3815
3816               x12=c(1,i)-c(1,j)
3817               y12=c(2,i)-c(2,j)
3818               z12=c(3,i)-c(3,j)
3819               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3820 c              write (iout,*) k,i,j,distal,dist2_cut
3821
3822             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3823      &            .and. distal.le.dist2_cut ) then
3824
3825               ii=ii+1
3826               ii_in_use(ii)=1
3827               l_homo(k,ii)=.true.
3828
3829 c             write (iout,*) "k",k
3830 c             write (iout,*) "i",i," j",j," constr_homology",
3831 c     &                       constr_homology
3832               ires_homo(ii)=i+chain_border1(1,igr)-1
3833               jres_homo(ii)=j+chain_border1(1,igr)-1
3834               odl(k,ii)=distal
3835               if (read2sigma) then
3836                 sigma_odl(k,ii)=0
3837                 do ik=i,j
3838                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3839                 enddo
3840                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3841                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3842      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3843               else
3844                 if (odl(k,ii).le.dist_cut) then
3845                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3846                 else
3847 #ifdef OLDSIGMA
3848                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3849      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3850 #else
3851                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3852      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3853 #endif
3854                 endif
3855               endif
3856               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3857             else
3858               ii=ii+1
3859 c              l_homo(k,ii)=.false.
3860             endif
3861             enddo
3862           enddo
3863         lim_odl=ii
3864         endif
3865 c
3866 c     Theta, dihedral and SC retraints
3867 c
3868         if (waga_angle.gt.0.0d0) then
3869           do i = nnt_chain+3,nct_chain
3870             iii=i+chain_border1(1,igr)-1
3871             if (idomain(k,i).eq.0) then
3872 c               sigma_dih(k,i)=0.0
3873                cycle
3874             endif
3875             dih(k,iii)=phiref(i)
3876             sigma_dih(k,iii)=
3877      &         (rescore(k,i)+rescore(k,i-1)+
3878      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3879 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3880 c     &       " sigma_dihed",sigma_dih(k,i)
3881             if (sigma_dih(k,iii).ne.0)
3882      &       sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
3883           enddo
3884 c          lim_dih=nct-nnt-2
3885         endif
3886
3887         if (waga_theta.gt.0.0d0) then
3888           do i = nnt_chain+2,nct_chain
3889              iii=i+chain_border1(1,igr)-1
3890              if (idomain(k,i).eq.0) then
3891 c              sigma_theta(k,i)=0.0
3892               cycle
3893              endif
3894              thetatpl(k,iii)=thetaref(i)
3895              sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+
3896      &                        rescore(k,i-2))/3.0
3897              if (sigma_theta(k,iii).ne.0)
3898      &       sigma_theta(k,iii)=1.0d0/
3899      &       (sigma_theta(k,iii)*sigma_theta(k,iii))
3900           enddo
3901         endif
3902
3903         if (waga_d.gt.0.0d0) then
3904           do i = nnt_chain,nct_chain
3905              iii=i+chain_border1(1,igr)-1
3906                if (itype(i).eq.10) cycle
3907                if (idomain(k,i).eq.0 ) then
3908 c                  sigma_d(k,i)=0.0
3909                   cycle
3910                endif
3911                xxtpl(k,iii)=xxref(i)
3912                yytpl(k,iii)=yyref(i)
3913                zztpl(k,iii)=zzref(i)
3914                sigma_d(k,iii)=rescore(k,i)
3915                if (sigma_d(k,iii).ne.0)
3916      &          sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
3917 c               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3918           enddo
3919         endif
3920       enddo ! l
3921       enddo ! ll
3922 c
3923 c remove distance restraints not used in any model from the list
3924 c shift data in all arrays
3925 c
3926 c      write (iout,*) "ii_old",ii_old
3927       if (waga_dist.ne.0.0d0) then
3928 #ifdef DEBUG
3929        write (iout,*) "Distance restraints from templates"
3930        do iii=1,lim_odl
3931        write(iout,'(4i5,100(2f8.2,1x,l1,4x))')
3932      &  iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii),
3933      &  (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii),
3934      &  ki=1,constr_homology)
3935        enddo
3936 #endif
3937         ii=ii_old
3938         liiflag=.true.
3939         lfirst=.true.
3940         do i=nnt_chain,nct_chain-2
3941          do j=i+2,nct_chain
3942           ii=ii+1
3943 c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3944 c     &            .and. distal.le.dist2_cut ) then
3945 c          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
3946 c          call flush(iout)
3947           if (ii_in_use(ii).eq.0.and.liiflag.or.
3948      &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3949             liiflag=.false.
3950             i10=ii
3951             if (lfirst) then
3952               lfirst=.false.
3953               iistart=ii
3954             else
3955               if(i10.eq.lim_odl) i10=i10+1
3956               do ki=0,i10-i01-1
3957                ires_homo(iistart+ki)=ires_homo(ki+i01)
3958                jres_homo(iistart+ki)=jres_homo(ki+i01)
3959                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3960                do k=1,constr_homology
3961                 odl(k,iistart+ki)=odl(k,ki+i01)
3962                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3963                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3964                enddo
3965               enddo
3966               iistart=iistart+i10-i01
3967             endif
3968           endif
3969           if (ii_in_use(ii).ne.0.and..not.liiflag) then
3970              i01=ii
3971              liiflag=.true.
3972           endif
3973          enddo
3974         enddo
3975         lim_odl=iistart-1
3976       endif
3977
3978       lll=lim_odl-ii_old
3979
3980       do i=2,nequiv(igr)
3981
3982         ichain=iequiv(i,igr)
3983
3984         do j=nnt_chain,nct_chain
3985           jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
3986           do k=1,constr_homology
3987             dih(k,jj)=dih(k,j)
3988             sigma_dih(k,jj)=sigma_dih(k,j)
3989             thetatpl(k,jj)=thetatpl(k,j)
3990             sigma_theta(k,jj)=sigma_theta(k,j)
3991             xxtpl(k,jj)=xxtpl(k,j)
3992             yytpl(k,jj)=yytpl(k,j)
3993             zztpl(k,jj)=zztpl(k,j)
3994             sigma_d(k,jj)=sigma_d(k,j)
3995           enddo
3996         enddo
3997
3998         jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
3999 c        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
4000         do j=ii_old+1,lim_odl
4001           ires_homo(j+lll)=ires_homo(j)+jj
4002           jres_homo(j+lll)=jres_homo(j)+jj
4003           do k=1,constr_homology
4004             odl(k,j+lll)=odl(k,j)
4005             sigma_odl(k,j+lll)=sigma_odl(k,j)
4006             l_homo(k,j+lll)=l_homo(k,j)
4007           enddo
4008         enddo
4009
4010         ii_old=ii_old+lll
4011         lim_odl=lim_odl+lll
4012
4013       enddo
4014
4015       ENDDO ! IGR
4016
4017       if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
4018       nres=nres_temp
4019       itype=itype_temp
4020
4021       return
4022    10 stop "Error in fragment file"
4023       end