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