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 c      write (iout,*) "iranconf",iranconf," extconf",extconf,
1219 c     & " start_from_models",start_from_model
1220       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1221      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1222      &    modecalc.ne.10) then
1223 C If input structure hasn't been supplied from the PDB file read or generate
1224 C initial geometry.
1225         if (iranconf.eq.0 .and. .not. extconf .and. .not.
1226      &     start_from_model) then
1227           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1228      &     write (iout,'(a)') 'Initial geometry will be read in.'
1229           if (read_cart) then
1230             read(inp,'(8f10.5)',end=36,err=36)
1231      &       ((c(l,k),l=1,3),k=1,nres),
1232      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1233 c            write (iout,*) "Exit READ_CART"
1234 c            write (iout,'(8f10.5)')
1235 c     &       ((c(l,k),l=1,3),k=1,nres),
1236 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1237             call cartprint
1238             do i=1,nres-1
1239               do j=1,3
1240                 dc(j,i)=c(j,i+1)-c(j,i)
1241 c                dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1242               enddo
1243             enddo
1244             do i=nnt,nct
1245               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1246                 do j=1,3
1247                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1248 c                  dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1249                 enddo
1250               else
1251                 do j=1,3
1252                   dc(j,i+nres)=0.0d0
1253 c                  dc_norm(j,i+nres)=0.0d0
1254                 enddo
1255               endif
1256             enddo
1257             call int_from_cart1(.true.)
1258             write (iout,*) "Finish INT_TO_CART"
1259 c      write (iout,*) "DC"
1260 c      do i=1,nres
1261 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1262 c     &    (dc(j,i+nres),j=1,3)
1263 c      enddo
1264 c            call cartprint
1265 c            call setup_var
1266 c            return
1267           else
1268             call read_angles(inp,*36)
1269             call bond_regular
1270             call chainbuild_extconf
1271           endif
1272           goto 37
1273    36     write (iout,'(a)') 'Error reading angle file.'
1274 #ifdef MPI
1275           call mpi_finalize( MPI_COMM_WORLD,IERR )
1276 #endif
1277           stop 'Error reading angle file.'
1278    37     continue
1279         else if (extconf) then
1280           if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1281      &      write (iout,'(a)') 'Extended chain initial geometry.'
1282           do i=3,nres
1283             theta(i)=90d0*deg2rad
1284           enddo
1285           do i=4,nres
1286             phi(i)=180d0*deg2rad
1287           enddo
1288           do i=2,nres-1
1289             alph(i)=110d0*deg2rad
1290           enddo
1291           do i=2,nres-1
1292             omeg(i)=-120d0*deg2rad
1293             if (itype(i).le.0) omeg(i)=-omeg(i)
1294           enddo
1295           call bond_regular
1296           call chainbuild_extconf
1297         else
1298           if(me.eq.king.or..not.out1file)
1299      &     write (iout,'(a)') 'Random-generated initial geometry.'
1300           call bond_regular
1301 #ifdef MPI
1302           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1303      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1304 #endif
1305             do itrial=1,100
1306               itmp=1
1307               call gen_rand_conf(itmp,*30)
1308               goto 40
1309    30         write (iout,*) 'Failed to generate random conformation',
1310      &          ', itrial=',itrial
1311               write (*,*) 'Processor:',me,
1312      &          ' Failed to generate random conformation',
1313      &          ' itrial=',itrial
1314               call intout
1315
1316 #ifdef AIX
1317               call flush_(iout)
1318 #else
1319               call flush(iout)
1320 #endif
1321             enddo
1322             write (iout,'(a,i3,a)') 'Processor:',me,
1323      &        ' error in generating random conformation.'
1324             write (*,'(a,i3,a)') 'Processor:',me,
1325      &        ' error in generating random conformation.'
1326             call flush(iout)
1327 #ifdef MPI
1328             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1329    40       continue
1330           endif
1331 #else
1332           write (*,'(a)')
1333      &      ' error in generating random conformation.'
1334           stop
1335    40     continue
1336 #endif
1337         endif
1338       elseif (modecalc.eq.4) then
1339         read (inp,'(a)') intinname
1340         open (intin,file=intinname,status='old',err=333)
1341         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1342      &  write (iout,'(a)') 'intinname',intinname
1343         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1344         goto 334
1345   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1346 #ifdef MPI
1347         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1348 #endif
1349         stop 'Error opening angle file.'
1350   334   continue
1351
1352       endif
1353 C Generate distance constraints, if the PDB structure is to be regularized.
1354       if (nthread.gt.0) then
1355         call read_threadbase
1356       endif
1357       call setup_var
1358       if (me.eq.king .or. .not. out1file)
1359      & call intout
1360       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1361         write (iout,'(/a,i3,a)')
1362      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1363         write (iout,'(20i4)') (iss(i),i=1,ns)
1364        if (dyn_ss) then
1365           write(iout,*)"Running with dynamic disulfide-bond formation"
1366        else
1367         write (iout,'(/a/)') 'Pre-formed links are:'
1368         do i=1,nss
1369           i1=ihpb(i)-nres
1370           i2=jhpb(i)-nres
1371           it1=itype(i1)
1372           it2=itype(i2)
1373           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1374      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1375      &    ebr,forcon(i)
1376         enddo
1377         write (iout,'(a)')
1378        endif
1379       endif
1380       if (ns.gt.0.and.dyn_ss) then
1381           do i=nss+1,nhpb
1382             ihpb(i-nss)=ihpb(i)
1383             jhpb(i-nss)=jhpb(i)
1384             forcon(i-nss)=forcon(i)
1385             dhpb(i-nss)=dhpb(i)
1386           enddo
1387           nhpb=nhpb-nss
1388           nss=0
1389           call hpb_partition
1390           do i=1,ns
1391             dyn_ss_mask(iss(i))=.true.
1392           enddo
1393       endif
1394 c --- test
1395 c      call cartprint
1396 c      write (iout,*) "DC"
1397 c      do i=1,nres
1398 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1399 c     &    (dc(j,i+nres),j=1,3)
1400 c      enddo
1401       if (i2ndstr.gt.0) call secstrp2dihc
1402 c      call geom_to_var(nvar,x)
1403 c      call etotal(energia(0))
1404 c      call enerprint(energia(0))
1405 c      call briefout(0,etot)
1406 c      stop
1407 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1408 cd    write (iout,'(a)') 'Variable list:'
1409 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1410 #ifdef MPI
1411       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1412      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1413      &  'Processor',myrank,': end reading molecular data.'
1414 #endif
1415 c      print *,"A TU?"
1416       return
1417       end
1418 c--------------------------------------------------------------------------
1419       logical function seq_comp(itypea,itypeb,length)
1420       implicit none
1421       integer length,itypea(length),itypeb(length)
1422       integer i
1423       do i=1,length
1424         if (itypea(i).ne.itypeb(i)) then
1425           seq_comp=.false.
1426           return
1427         endif
1428       enddo
1429       seq_comp=.true.
1430       return
1431       end
1432 c-----------------------------------------------------------------------------
1433       subroutine read_bridge
1434 C Read information about disulfide bridges.
1435       implicit none
1436       include 'DIMENSIONS'
1437 #ifdef MPI
1438       include 'mpif.h'
1439       integer ierror
1440 #endif
1441       include 'COMMON.IOUNITS'
1442       include 'COMMON.GEO'
1443       include 'COMMON.VAR'
1444       include 'COMMON.INTERACT'
1445       include 'COMMON.LOCAL'
1446       include 'COMMON.NAMES'
1447       include 'COMMON.CHAIN'
1448       include 'COMMON.FFIELD'
1449       include 'COMMON.SBRIDGE'
1450       include 'COMMON.HEADER'
1451       include 'COMMON.CONTROL'
1452       include 'COMMON.DBASE'
1453       include 'COMMON.THREAD'
1454       include 'COMMON.TIME1'
1455       include 'COMMON.SETUP'
1456       integer i,j
1457 C Read bridging residues.
1458       read (inp,*) ns,(iss(i),i=1,ns)
1459 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1460 c     numbers
1461       icys=0
1462       do i=1,ns
1463         icys(iss(i))=i
1464       enddo
1465 c      print *,'ns=',ns
1466       if(me.eq.king.or..not.out1file)
1467      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1468 C Check whether the specified bridging residues are cystines.
1469       do i=1,ns
1470         if (itype(iss(i)).ne.1) then
1471           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1472      &   'Do you REALLY think that the residue ',
1473      &    restyp(itype(iss(i))),i,
1474      &   ' can form a disulfide bridge?!!!'
1475           write (*,'(2a,i3,a)')
1476      &   'Do you REALLY think that the residue ',
1477      &    restyp(itype(iss(i))),i,
1478      &   ' can form a disulfide bridge?!!!'
1479 #ifdef MPI
1480          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1481          stop
1482 #endif
1483         endif
1484       enddo
1485 C Read preformed bridges.
1486       if (ns.gt.0) then
1487       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1488       if(fg_rank.eq.0)
1489      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1490       if (nss.gt.0) then
1491         nhpb=nss
1492 C Check if the residues involved in bridges are in the specified list of
1493 C bridging residues.
1494         do i=1,nss
1495           do j=1,i-1
1496             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1497      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1498               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1499      &      ' contains residues present in other pairs.'
1500               write (*,'(a,i3,a)') 'Disulfide pair',i,
1501      &      ' contains residues present in other pairs.'
1502 #ifdef MPI
1503               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1504               stop
1505 #endif
1506             endif
1507           enddo
1508           do j=1,ns
1509             if (ihpb(i).eq.iss(j)) goto 10
1510           enddo
1511           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1512    10     continue
1513           do j=1,ns
1514             if (jhpb(i).eq.iss(j)) goto 20
1515           enddo
1516           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1517    20     continue
1518 c          dhpb(i)=dbr
1519 c          forcon(i)=fbr
1520         enddo
1521         do i=1,nss
1522           ihpb(i)=ihpb(i)+nres
1523           jhpb(i)=jhpb(i)+nres
1524         enddo
1525       endif
1526       endif
1527       return
1528       end
1529 c----------------------------------------------------------------------------
1530       subroutine read_x(kanal,*)
1531       implicit none
1532       include 'DIMENSIONS'
1533       include 'COMMON.GEO'
1534       include 'COMMON.VAR'
1535       include 'COMMON.CHAIN'
1536       include 'COMMON.IOUNITS'
1537       include 'COMMON.CONTROL'
1538       include 'COMMON.LOCAL'
1539       include 'COMMON.INTERACT'
1540       integer i,j,k,l,kanal
1541 c Read coordinates from input
1542 c
1543       read(kanal,'(8f10.5)',end=10,err=10)
1544      &  ((c(l,k),l=1,3),k=1,nres),
1545      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1546       do j=1,3
1547         c(j,nres+1)=c(j,1)
1548         c(j,2*nres)=c(j,nres)
1549       enddo
1550       call int_from_cart1(.false.)
1551       do i=1,nres-1
1552         do j=1,3
1553           dc(j,i)=c(j,i+1)-c(j,i)
1554           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1555         enddo
1556       enddo
1557       do i=nnt,nct
1558         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1559           do j=1,3
1560             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1561             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1562           enddo
1563         endif
1564       enddo
1565
1566       return
1567    10 return1
1568       end
1569 c----------------------------------------------------------------------------
1570       subroutine read_threadbase
1571       implicit none
1572       include 'DIMENSIONS'
1573       include 'COMMON.IOUNITS'
1574       include 'COMMON.GEO'
1575       include 'COMMON.VAR'
1576       include 'COMMON.INTERACT'
1577       include 'COMMON.LOCAL'
1578       include 'COMMON.NAMES'
1579       include 'COMMON.CHAIN'
1580       include 'COMMON.FFIELD'
1581       include 'COMMON.SBRIDGE'
1582       include 'COMMON.HEADER'
1583       include 'COMMON.CONTROL'
1584       include 'COMMON.DBASE'
1585       include 'COMMON.THREAD'
1586       include 'COMMON.TIME1'
1587       integer i,j,k
1588       double precision dist
1589 C Read pattern database for threading.
1590       read (icbase,*) nseq
1591       do i=1,nseq
1592         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1593      &   nres_base(2,i),nres_base(3,i)
1594         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1595      &   nres_base(1,i))
1596 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1597 c    &   nres_base(2,i),nres_base(3,i)
1598 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1599 c    &   nres_base(1,i))
1600       enddo
1601       close (icbase)
1602       if (weidis.eq.0.0D0) weidis=0.1D0
1603       do i=nnt,nct
1604         do j=i+2,nct
1605           nhpb=nhpb+1
1606           ihpb(nhpb)=i
1607           jhpb(nhpb)=j
1608           forcon(nhpb)=weidis
1609         enddo
1610       enddo
1611       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1612       write (iout,'(a,i5)') 'nexcl: ',nexcl
1613       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1614       return
1615       end
1616 c------------------------------------------------------------------------------
1617       subroutine setup_var
1618       implicit none
1619       include 'DIMENSIONS'
1620       include 'COMMON.IOUNITS'
1621       include 'COMMON.GEO'
1622       include 'COMMON.VAR'
1623       include 'COMMON.INTERACT'
1624       include 'COMMON.LOCAL'
1625       include 'COMMON.NAMES'
1626       include 'COMMON.CHAIN'
1627       include 'COMMON.FFIELD'
1628       include 'COMMON.SBRIDGE'
1629       include 'COMMON.HEADER'
1630       include 'COMMON.CONTROL'
1631       include 'COMMON.DBASE'
1632       include 'COMMON.THREAD'
1633       include 'COMMON.TIME1'
1634       integer i
1635 C Set up variable list.
1636       ntheta=nres-2
1637       nphi=nres-3
1638       nvar=ntheta+nphi
1639       nside=0
1640       do i=2,nres-1
1641         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1642           nside=nside+1
1643           ialph(i,1)=nvar+nside
1644           ialph(nside,2)=i
1645         endif
1646       enddo
1647       if (indphi.gt.0) then
1648         nvar=nphi
1649       else if (indback.gt.0) then
1650         nvar=nphi+ntheta
1651       else
1652         nvar=nvar+2*nside
1653       endif
1654       return
1655       end
1656 c----------------------------------------------------------------------------
1657       subroutine gen_dist_constr
1658 C Generate CA distance constraints.
1659       implicit none
1660       include 'DIMENSIONS'
1661       include 'COMMON.IOUNITS'
1662       include 'COMMON.GEO'
1663       include 'COMMON.VAR'
1664       include 'COMMON.INTERACT'
1665       include 'COMMON.LOCAL'
1666       include 'COMMON.NAMES'
1667       include 'COMMON.CHAIN'
1668       include 'COMMON.FFIELD'
1669       include 'COMMON.SBRIDGE'
1670       include 'COMMON.HEADER'
1671       include 'COMMON.CONTROL'
1672       include 'COMMON.DBASE'
1673       include 'COMMON.THREAD'
1674       include 'COMMON.SPLITELE'
1675       include 'COMMON.TIME1'
1676       integer i,j,itype_pdb(maxres)
1677       common /pizda/ itype_pdb
1678       double precision dd
1679       double precision dist
1680       character*2 iden
1681 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1682 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1683 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1684 cd     & ' nsup',nsup
1685       do i=nstart_sup,nstart_sup+nsup-1
1686 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1687 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1688         do j=i+2,nstart_sup+nsup-1
1689 c 5/24/2020 Adam: Cutoff included to reduce array size
1690           dd = dist(i,j)
1691           if (dd.gt.r_cut_int) cycle
1692           nhpb=nhpb+1
1693           ihpb(nhpb)=i+nstart_seq-nstart_sup
1694           jhpb(nhpb)=j+nstart_seq-nstart_sup
1695           forcon(nhpb)=weidis
1696           dhpb(nhpb)=dd
1697         enddo
1698       enddo
1699 cd      write (iout,'(a)') 'Distance constraints:'
1700 cd      do i=nss+1,nhpb
1701 cd        ii=ihpb(i)
1702 cd        jj=jhpb(i)
1703 cd        iden='CA'
1704 cd        if (ii.gt.nres) then
1705 cd          iden='SC'
1706 cd          ii=ii-nres
1707 cd          jj=jj-nres
1708 cd        endif
1709 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1710 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1711 cd     &  dhpb(i),forcon(i)
1712 cd      enddo
1713       return
1714       end
1715 c----------------------------------------------------------------------------
1716       subroutine map_read
1717       implicit none
1718       include 'DIMENSIONS'
1719       include 'COMMON.MAP'
1720       include 'COMMON.IOUNITS'
1721       integer imap
1722       character*3 angid(4) /'THE','PHI','ALP','OME'/
1723       character*80 mapcard,ucase
1724       do imap=1,nmap
1725         read (inp,'(a)') mapcard
1726         mapcard=ucase(mapcard)
1727         if (index(mapcard,'PHI').gt.0) then
1728           kang(imap)=1
1729         else if (index(mapcard,'THE').gt.0) then
1730           kang(imap)=2
1731         else if (index(mapcard,'ALP').gt.0) then
1732           kang(imap)=3
1733         else if (index(mapcard,'OME').gt.0) then
1734           kang(imap)=4
1735         else
1736           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1737           stop 'Error - illegal variable spec in MAP card.'
1738         endif
1739         call readi (mapcard,'RES1',res1(imap),0)
1740         call readi (mapcard,'RES2',res2(imap),0)
1741         if (res1(imap).eq.0) then
1742           res1(imap)=res2(imap)
1743         else if (res2(imap).eq.0) then
1744           res2(imap)=res1(imap)
1745         endif
1746         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1747           write (iout,'(a)')
1748      &    'Error - illegal definition of variable group in MAP.'
1749           stop 'Error - illegal definition of variable group in MAP.'
1750         endif
1751         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1752         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1753         call readi(mapcard,'NSTEP',nstep(imap),0)
1754         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1755           write (iout,'(a)')
1756      &     'Illegal boundary and/or step size specification in MAP.'
1757           stop 'Illegal boundary and/or step size specification in MAP.'
1758         endif
1759       enddo ! imap
1760       return
1761       end
1762 c----------------------------------------------------------------------------
1763       subroutine csaread
1764       implicit none
1765       include 'DIMENSIONS'
1766       include 'COMMON.IOUNITS'
1767       include 'COMMON.GEO'
1768       include 'COMMON.CSA'
1769       include 'COMMON.BANK'
1770       include 'COMMON.CONTROL'
1771       character*80 ucase
1772       character*620 mcmcard
1773       call card_concat(mcmcard)
1774
1775       call readi(mcmcard,'NCONF',nconf,50)
1776       call readi(mcmcard,'NADD',nadd,0)
1777       call readi(mcmcard,'JSTART',jstart,1)
1778       call readi(mcmcard,'JEND',jend,1)
1779       call readi(mcmcard,'NSTMAX',nstmax,500000)
1780       call readi(mcmcard,'N0',n0,1)
1781       call readi(mcmcard,'N1',n1,6)
1782       call readi(mcmcard,'N2',n2,4)
1783       call readi(mcmcard,'N3',n3,0)
1784       call readi(mcmcard,'N4',n4,0)
1785       call readi(mcmcard,'N5',n5,0)
1786       call readi(mcmcard,'N6',n6,10)
1787       call readi(mcmcard,'N7',n7,0)
1788       call readi(mcmcard,'N8',n8,0)
1789       call readi(mcmcard,'N9',n9,0)
1790       call readi(mcmcard,'N14',n14,0)
1791       call readi(mcmcard,'N15',n15,0)
1792       call readi(mcmcard,'N16',n16,0)
1793       call readi(mcmcard,'N17',n17,0)
1794       call readi(mcmcard,'N18',n18,0)
1795
1796       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1797
1798       call readi(mcmcard,'NDIFF',ndiff,2)
1799       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1800       call readi(mcmcard,'IS1',is1,1)
1801       call readi(mcmcard,'IS2',is2,8)
1802       call readi(mcmcard,'NRAN0',nran0,4)
1803       call readi(mcmcard,'NRAN1',nran1,2)
1804       call readi(mcmcard,'IRR',irr,1)
1805       call readi(mcmcard,'NSEED',nseed,20)
1806       call readi(mcmcard,'NTOTAL',ntotal,10000)
1807       call reada(mcmcard,'CUT1',cut1,2.0d0)
1808       call reada(mcmcard,'CUT2',cut2,5.0d0)
1809       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1810       call readi(mcmcard,'ICMAX',icmax,3)
1811       call readi(mcmcard,'IRESTART',irestart,0)
1812 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1813       ntbankm=0
1814 c!bankt
1815       call reada(mcmcard,'DELE',dele,20.0d0)
1816       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1817       call readi(mcmcard,'IREF',iref,0)
1818       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1819       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1820       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1821       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1822       write (iout,*) "NCONF_IN",nconf_in
1823       return
1824       end
1825 c----------------------------------------------------------------------------
1826       subroutine mcmread
1827       implicit none
1828       include 'DIMENSIONS'
1829       include 'COMMON.MCM'
1830       include 'COMMON.MCE'
1831       include 'COMMON.IOUNITS'
1832       character*80 ucase
1833       character*320 mcmcard
1834       integer i
1835       call card_concat(mcmcard)
1836       call readi(mcmcard,'MAXACC',maxacc,100)
1837       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1838       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1839       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1840       call readi(mcmcard,'MAXREPM',maxrepm,200)
1841       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1842       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1843       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1844       call reada(mcmcard,'E_UP',e_up,5.0D0)
1845       call reada(mcmcard,'DELTE',delte,0.1D0)
1846       call readi(mcmcard,'NSWEEP',nsweep,5)
1847       call readi(mcmcard,'NSTEPH',nsteph,0)
1848       call readi(mcmcard,'NSTEPC',nstepc,0)
1849       call reada(mcmcard,'TMIN',tmin,298.0D0)
1850       call reada(mcmcard,'TMAX',tmax,298.0D0)
1851       call readi(mcmcard,'NWINDOW',nwindow,0)
1852       call readi(mcmcard,'PRINT_MC',print_mc,0)
1853       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1854       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1855       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1856       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1857       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1858       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1859       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1860       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1861       if (nwindow.gt.0) then
1862         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1863         do i=1,nwindow
1864           winlen(i)=winend(i)-winstart(i)+1
1865         enddo
1866       endif
1867       if (tmax.lt.tmin) tmax=tmin
1868       if (tmax.eq.tmin) then
1869         nstepc=0
1870         nsteph=0
1871       endif
1872       if (nstepc.gt.0 .and. nsteph.gt.0) then
1873         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1874         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1875       endif
1876 C Probabilities of different move types
1877       sumpro_type(0)=0.0D0
1878       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1879       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1880       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1881       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1882       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1883       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1884       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1885       do i=1,MaxMoveType
1886         print *,'i',i,' sumprotype',sumpro_type(i)
1887         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1888         print *,'i',i,' sumprotype',sumpro_type(i)
1889       enddo
1890       return
1891       end
1892 c----------------------------------------------------------------------------
1893       subroutine read_minim
1894       implicit none
1895       include 'DIMENSIONS'
1896       include 'COMMON.MINIM'
1897       include 'COMMON.IOUNITS'
1898       character*80 ucase
1899       character*320 minimcard
1900       call card_concat(minimcard)
1901       call readi(minimcard,'MAXMIN',maxmin,2000)
1902       call readi(minimcard,'MAXFUN',maxfun,5000)
1903       call readi(minimcard,'MINMIN',minmin,maxmin)
1904       call readi(minimcard,'MINFUN',minfun,maxmin)
1905       call reada(minimcard,'TOLF',tolf,1.0D-2)
1906       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1907       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1908       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1909       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1910       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1911      &         'Options in energy minimization:'
1912       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1913      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1914      & 'MinMin:',MinMin,' MinFun:',MinFun,
1915      & ' TolF:',TolF,' RTolF:',RTolF
1916       return
1917       end
1918 c----------------------------------------------------------------------------
1919       subroutine read_angles(kanal,*)
1920       implicit none
1921       include 'DIMENSIONS'
1922       include 'COMMON.GEO'
1923       include 'COMMON.VAR'
1924       include 'COMMON.CHAIN'
1925       include 'COMMON.IOUNITS'
1926       include 'COMMON.CONTROL'
1927       integer i,kanal
1928 c Read angles from input
1929 c
1930        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1931        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1932        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1933        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1934
1935        do i=1,nres
1936 c 9/7/01 avoid 180 deg valence angle
1937         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1938 c
1939         theta(i)=deg2rad*theta(i)
1940         phi(i)=deg2rad*phi(i)
1941         alph(i)=deg2rad*alph(i)
1942         omeg(i)=deg2rad*omeg(i)
1943        enddo
1944       return
1945    10 return1
1946       end
1947 c----------------------------------------------------------------------------
1948       subroutine reada(rekord,lancuch,wartosc,default)
1949       implicit none
1950       character*(*) rekord,lancuch
1951       double precision wartosc,default
1952       integer ilen,iread
1953       external ilen
1954       iread=index(rekord,lancuch)
1955       if (iread.eq.0) then
1956         wartosc=default
1957         return
1958       endif
1959       iread=iread+ilen(lancuch)+1
1960       read (rekord(iread:),*,err=10,end=10) wartosc
1961       return
1962   10  wartosc=default
1963       return
1964       end
1965 c----------------------------------------------------------------------------
1966       subroutine readi(rekord,lancuch,wartosc,default)
1967       implicit none
1968       character*(*) rekord,lancuch
1969       integer wartosc,default
1970       integer ilen,iread
1971       external ilen
1972       iread=index(rekord,lancuch)
1973       if (iread.eq.0) then
1974         wartosc=default
1975         return
1976       endif
1977       iread=iread+ilen(lancuch)+1
1978       read (rekord(iread:),*,err=10,end=10) wartosc
1979       return
1980   10  wartosc=default
1981       return
1982       end
1983 c----------------------------------------------------------------------------
1984       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1985       implicit none
1986       integer dim,i
1987       integer tablica(dim),default
1988       character*(*) rekord,lancuch
1989       character*80 aux
1990       integer ilen,iread
1991       external ilen
1992       do i=1,dim
1993         tablica(i)=default
1994       enddo
1995       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1996       if (iread.eq.0) return
1997       iread=iread+ilen(lancuch)+1
1998       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1999    10 return
2000       end
2001 c----------------------------------------------------------------------------
2002       subroutine multreada(rekord,lancuch,tablica,dim,default)
2003       implicit none
2004       integer dim,i
2005       double precision tablica(dim),default
2006       character*(*) rekord,lancuch
2007       character*80 aux
2008       integer ilen,iread
2009       external ilen
2010       do i=1,dim
2011         tablica(i)=default
2012       enddo
2013       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2014       if (iread.eq.0) return
2015       iread=iread+ilen(lancuch)+1
2016       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2017    10 return
2018       end
2019 c----------------------------------------------------------------------------
2020       subroutine openunits
2021       implicit none
2022       include 'DIMENSIONS'
2023 #ifdef MPI
2024       include 'mpif.h'
2025       integer ierror
2026       character*16 form,nodename
2027       integer nodelen
2028 #endif
2029       include 'COMMON.SETUP'
2030       include 'COMMON.IOUNITS'
2031       include 'COMMON.MD'
2032       include 'COMMON.CONTROL'
2033       integer lenpre,lenpot,ilen,lentmp,npos
2034       external ilen
2035       character*3 out1file_text,ucase
2036       character*3 ll
2037       external ucase
2038       tmpdir=""
2039       out1file_text="YES"
2040 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2041       call getenv_loc("PREFIX",prefix)
2042       pref_orig = prefix
2043       call getenv_loc("POT",pot)
2044       call getenv_loc("DIRTMP",tmpdir)
2045       call getenv_loc("CURDIR",curdir)
2046       call getenv_loc("OUT1FILE",out1file_text)
2047 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2048       out1file_text=ucase(out1file_text)
2049       if (out1file_text(1:1).eq."Y") then
2050         out1file=.true.
2051       else
2052         out1file=fg_rank.gt.0
2053       endif
2054       lenpre=ilen(prefix)
2055       lenpot=ilen(pot)
2056       lentmp=ilen(tmpdir)
2057       if (lentmp.gt.0) then
2058           write (*,'(80(1h!))')
2059           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2060           write (*,'(80(1h!))')
2061           write (*,*)"All output files will be on node /tmp directory."
2062 #ifdef MPI
2063         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2064         if (me.eq.king) then
2065           write (*,*) "The master node is ",nodename
2066         else if (fg_rank.eq.0) then
2067           write (*,*) "I am the CG slave node ",nodename
2068         else
2069           write (*,*) "I am the FG slave node ",nodename
2070         endif
2071 #endif
2072         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2073         lenpre = lentmp+lenpre+1
2074       endif
2075       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2076 C Get the names and open the input files
2077 #if defined(WINIFL) || defined(WINPGI)
2078       open(1,file=pref_orig(:ilen(pref_orig))//
2079      &  '.inp',status='old',readonly,shared)
2080        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2081 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2082 C Get parameter filenames and open the parameter files.
2083       call getenv_loc('BONDPAR',bondname)
2084       open (ibond,file=bondname,status='old',readonly,shared)
2085       call getenv_loc('THETPAR',thetname)
2086       open (ithep,file=thetname,status='old',readonly,shared)
2087       call getenv_loc('ROTPAR',rotname)
2088       open (irotam,file=rotname,status='old',readonly,shared)
2089       call getenv_loc('TORPAR',torname)
2090       open (itorp,file=torname,status='old',readonly,shared)
2091       call getenv_loc('TORDPAR',tordname)
2092       open (itordp,file=tordname,status='old',readonly,shared)
2093       call getenv_loc('FOURIER',fouriername)
2094       open (ifourier,file=fouriername,status='old',readonly,shared)
2095       call getenv_loc('ELEPAR',elename)
2096       open (ielep,file=elename,status='old',readonly,shared)
2097       call getenv_loc('SIDEPAR',sidename)
2098       open (isidep,file=sidename,status='old',readonly,shared)
2099       call getenv_loc('LIPTRANPAR',liptranname)
2100       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2101 #elif (defined CRAY) || (defined AIX)
2102       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2103      &  action='read')
2104 c      print *,"Processor",myrank," opened file 1"
2105       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2106 c      print *,"Processor",myrank," opened file 9"
2107 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2108 C Get parameter filenames and open the parameter files.
2109       call getenv_loc('BONDPAR',bondname)
2110       open (ibond,file=bondname,status='old',action='read')
2111 c      print *,"Processor",myrank," opened file IBOND"
2112       call getenv_loc('THETPAR',thetname)
2113       open (ithep,file=thetname,status='old',action='read')
2114 #ifndef CRYST_THETA
2115       call getenv_loc('THETPARPDB',thetname_pdb)
2116       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2117 #endif
2118 c      print *,"Processor",myrank," opened file ITHEP"
2119       call getenv_loc('ROTPAR',rotname)
2120       open (irotam,file=rotname,status='old',action='read')
2121 #ifndef CRYST_SC
2122       call getenv_loc('ROTPARPDB',rotname_pdb)
2123       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2124 #endif
2125 c      print *,"Processor",myrank," opened file IROTAM"
2126       call getenv_loc('TORPAR',torname)
2127       open (itorp,file=torname,status='old',action='read')
2128 c      print *,"Processor",myrank," opened file ITORP"
2129       call getenv_loc('TORDPAR',tordname)
2130       open (itordp,file=tordname,status='old',action='read')
2131 c      print *,"Processor",myrank," opened file ITORDP"
2132       call getenv_loc('SCCORPAR',sccorname)
2133       open (isccor,file=sccorname,status='old',action='read')
2134 c      print *,"Processor",myrank," opened file ISCCOR"
2135       call getenv_loc('FOURIER',fouriername)
2136       open (ifourier,file=fouriername,status='old',action='read')
2137 c      print *,"Processor",myrank," opened file IFOURIER"
2138       call getenv_loc('ELEPAR',elename)
2139       open (ielep,file=elename,status='old',action='read')
2140 c      print *,"Processor",myrank," opened file IELEP"
2141       call getenv_loc('SIDEPAR',sidename)
2142       open (isidep,file=sidename,status='old',action='read')
2143       call getenv_loc('LIPTRANPAR',liptranname)
2144       open (iliptranpar,file=liptranname,status='old',action='read')
2145 c      print *,"Processor",myrank," opened file ISIDEP"
2146 c      print *,"Processor",myrank," opened parameter files"
2147 #elif (defined G77)
2148       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2149       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2150 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2151 C Get parameter filenames and open the parameter files.
2152       call getenv_loc('BONDPAR',bondname)
2153       open (ibond,file=bondname,status='old')
2154       call getenv_loc('THETPAR',thetname)
2155       open (ithep,file=thetname,status='old')
2156 #ifndef CRYST_THETA
2157       call getenv_loc('THETPARPDB',thetname_pdb)
2158       open (ithep_pdb,file=thetname_pdb,status='old')
2159 #endif
2160       call getenv_loc('ROTPAR',rotname)
2161       open (irotam,file=rotname,status='old')
2162 #ifndef CRYST_SC
2163       call getenv_loc('ROTPARPDB',rotname_pdb)
2164       open (irotam_pdb,file=rotname_pdb,status='old')
2165 #endif
2166       call getenv_loc('TORPAR',torname)
2167       open (itorp,file=torname,status='old')
2168       call getenv_loc('TORDPAR',tordname)
2169       open (itordp,file=tordname,status='old')
2170       call getenv_loc('SCCORPAR',sccorname)
2171       open (isccor,file=sccorname,status='old')
2172       call getenv_loc('FOURIER',fouriername)
2173       open (ifourier,file=fouriername,status='old')
2174       call getenv_loc('ELEPAR',elename)
2175       open (ielep,file=elename,status='old')
2176       call getenv_loc('SIDEPAR',sidename)
2177       open (isidep,file=sidename,status='old')
2178       call getenv_loc('LIPTRANPAR',liptranname)
2179       open (iliptranpar,file=liptranname,status='old')
2180 #else
2181       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2182      &  readonly)
2183        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2184 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2185 C Get parameter filenames and open the parameter files.
2186       call getenv_loc('BONDPAR',bondname)
2187       open (ibond,file=bondname,status='old',readonly)
2188       call getenv_loc('THETPAR',thetname)
2189       open (ithep,file=thetname,status='old',readonly)
2190       call getenv_loc('ROTPAR',rotname)
2191       open (irotam,file=rotname,status='old',readonly)
2192       call getenv_loc('TORPAR',torname)
2193       open (itorp,file=torname,status='old',readonly)
2194       call getenv_loc('TORDPAR',tordname)
2195       open (itordp,file=tordname,status='old',readonly)
2196       call getenv_loc('SCCORPAR',sccorname)
2197       open (isccor,file=sccorname,status='old',readonly)
2198 #ifndef CRYST_THETA
2199       call getenv_loc('THETPARPDB',thetname_pdb)
2200       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2201 #endif
2202       call getenv_loc('FOURIER',fouriername)
2203       open (ifourier,file=fouriername,status='old',readonly)
2204       call getenv_loc('ELEPAR',elename)
2205       open (ielep,file=elename,status='old',readonly)
2206       call getenv_loc('SIDEPAR',sidename)
2207       open (isidep,file=sidename,status='old',readonly)
2208       call getenv_loc('LIPTRANPAR',liptranname)
2209       open (iliptranpar,file=liptranname,status='old',action='read')
2210 #ifndef CRYST_SC
2211       call getenv_loc('ROTPARPDB',rotname_pdb)
2212       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2213 #endif
2214 #endif
2215 #ifdef TUBE
2216       call getenv_loc('TUBEPAR',tubename)
2217 #if defined(WINIFL) || defined(WINPGI)
2218       open (itube,file=tubename,status='old',readonly,shared)
2219 #elif (defined CRAY)  || (defined AIX)
2220       open (itube,file=tubename,status='old',action='read')
2221 #elif (defined G77)
2222       open (itube,file=tubename,status='old')
2223 #else
2224       open (itube,file=tubename,status='old',readonly)
2225 #endif
2226 #endif
2227 #ifndef OLDSCP
2228 C
2229 C 8/9/01 In the newest version SCp interaction constants are read from a file
2230 C Use -DOLDSCP to use hard-coded constants instead.
2231 C
2232       call getenv_loc('SCPPAR',scpname)
2233 #if defined(WINIFL) || defined(WINPGI)
2234       open (iscpp,file=scpname,status='old',readonly,shared)
2235 #elif (defined CRAY)  || (defined AIX)
2236       open (iscpp,file=scpname,status='old',action='read')
2237 #elif (defined G77)
2238       open (iscpp,file=scpname,status='old')
2239 #else
2240       open (iscpp,file=scpname,status='old',readonly)
2241 #endif
2242 #endif
2243       call getenv_loc('PATTERN',patname)
2244 #if defined(WINIFL) || defined(WINPGI)
2245       open (icbase,file=patname,status='old',readonly,shared)
2246 #elif (defined CRAY)  || (defined AIX)
2247       open (icbase,file=patname,status='old',action='read')
2248 #elif (defined G77)
2249       open (icbase,file=patname,status='old')
2250 #else
2251       open (icbase,file=patname,status='old',readonly)
2252 #endif
2253 #ifdef MPI
2254 C Open output file only for CG processes
2255 c      print *,"Processor",myrank," fg_rank",fg_rank
2256       if (fg_rank.eq.0) then
2257
2258       if (nodes.eq.1) then
2259         npos=3
2260       else
2261         npos = dlog10(dfloat(nodes-1))+1
2262       endif
2263       if (npos.lt.3) npos=3
2264       write (liczba,'(i1)') npos
2265       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2266      &  //')'
2267       write (liczba,form) me
2268       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2269      &  liczba(:ilen(liczba))
2270       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2271      &  //'.int'
2272       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2273      &  //'.pdb'
2274       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2275      &  liczba(:ilen(liczba))//'.mol2'
2276       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2277      &  liczba(:ilen(liczba))//'.stat'
2278       if (lentmp.gt.0)
2279      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2280      &      //liczba(:ilen(liczba))//'.stat')
2281       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2282      &  //'.rst'
2283       if(usampl) then
2284           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2285      & liczba(:ilen(liczba))//'.const'
2286       endif
2287
2288       endif
2289 #else
2290       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2291       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2292       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2293       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2294       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2295       if (lentmp.gt.0)
2296      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2297      &    //'.stat')
2298       rest2name=prefix(:ilen(prefix))//'.rst'
2299       if(usampl) then
2300          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2301       endif
2302 #endif
2303 #if defined(AIX) || defined(PGI) || defined(CRAY)
2304       if (me.eq.king .or. .not. out1file) then
2305          open(iout,file=outname,status='unknown')
2306       else
2307          open(iout,file="/dev/null",status="unknown")
2308       endif
2309 c#define DEBUG
2310 #ifdef DEBUG
2311       if (fg_rank.gt.0) then
2312         write (liczba,'(i3.3)') myrank/nfgtasks
2313         write (ll,'(bz,i3.3)') fg_rank
2314         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2315      &   status='unknown')
2316       endif
2317 #endif
2318 c#undef DEBUG
2319       if(me.eq.king) then
2320        open(igeom,file=intname,status='unknown',position='append')
2321        open(ipdb,file=pdbname,status='unknown')
2322        open(imol2,file=mol2name,status='unknown')
2323        open(istat,file=statname,status='unknown',position='append')
2324       else
2325 c1out       open(iout,file=outname,status='unknown')
2326       endif
2327 #else
2328       if (me.eq.king .or. .not.out1file)
2329      &    open(iout,file=outname,status='unknown')
2330 #ifdef DEBUG
2331       if (fg_rank.gt.0) then
2332         write (liczba,'(i3.3)') myrank/nfgtasks
2333         write (ll,'(bz,i3.3)') fg_rank
2334         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2335      &   status='unknown')
2336       endif
2337 #endif
2338       if(me.eq.king) then
2339        open(igeom,file=intname,status='unknown',access='append')
2340        open(ipdb,file=pdbname,status='unknown')
2341        open(imol2,file=mol2name,status='unknown')
2342        open(istat,file=statname,status='unknown',access='append')
2343       else
2344 c1out       open(iout,file=outname,status='unknown')
2345       endif
2346 #endif
2347       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2348       csa_seed=prefix(:lenpre)//'.CSA.seed'
2349       csa_history=prefix(:lenpre)//'.CSA.history'
2350       csa_bank=prefix(:lenpre)//'.CSA.bank'
2351       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2352       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2353       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2354 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2355       csa_int=prefix(:lenpre)//'.int'
2356       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2357       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2358       csa_in=prefix(:lenpre)//'.CSA.in'
2359 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2360 C Write file names
2361       if (me.eq.king)then
2362       write (iout,'(80(1h-))')
2363       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2364       write (iout,'(80(1h-))')
2365       write (iout,*) "Input file                      : ",
2366      &  pref_orig(:ilen(pref_orig))//'.inp'
2367       write (iout,*) "Output file                     : ",
2368      &  outname(:ilen(outname))
2369       write (iout,*)
2370       write (iout,*) "Sidechain potential file        : ",
2371      &  sidename(:ilen(sidename))
2372 #ifndef OLDSCP
2373       write (iout,*) "SCp potential file              : ",
2374      &  scpname(:ilen(scpname))
2375 #endif
2376       write (iout,*) "Electrostatic potential file    : ",
2377      &  elename(:ilen(elename))
2378       write (iout,*) "Cumulant coefficient file       : ",
2379      &  fouriername(:ilen(fouriername))
2380       write (iout,*) "Torsional parameter file        : ",
2381      &  torname(:ilen(torname))
2382       write (iout,*) "Double torsional parameter file : ",
2383      &  tordname(:ilen(tordname))
2384       write (iout,*) "SCCOR parameter file : ",
2385      &  sccorname(:ilen(sccorname))
2386       write (iout,*) "Bond & inertia constant file    : ",
2387      &  bondname(:ilen(bondname))
2388       write (iout,*) "Bending-torsion parameter file  : ",
2389      &  thetname(:ilen(thetname))
2390       write (iout,*) "Rotamer parameter file          : ",
2391      &  rotname(:ilen(rotname))
2392       write (iout,*) "Thetpdb parameter file          : ",
2393      &  thetname_pdb(:ilen(thetname_pdb))
2394       write (iout,*) "Threading database              : ",
2395      &  patname(:ilen(patname))
2396       if (lentmp.ne.0)
2397      &write (iout,*)" DIRTMP                          : ",
2398      &  tmpdir(:lentmp)
2399       write (iout,'(80(1h-))')
2400       endif
2401       return
2402       end
2403 c----------------------------------------------------------------------------
2404       subroutine card_concat(card)
2405       implicit real*8 (a-h,o-z)
2406       include 'DIMENSIONS'
2407       include 'COMMON.IOUNITS'
2408       character*(*) card
2409       character*80 karta,ucase
2410       external ilen
2411       read (inp,'(a)') karta
2412       karta=ucase(karta)
2413       card=' '
2414       do while (karta(80:80).eq.'&')
2415         card=card(:ilen(card)+1)//karta(:79)
2416         read (inp,'(a)') karta
2417         karta=ucase(karta)
2418       enddo
2419       card=card(:ilen(card)+1)//karta
2420       return
2421       end
2422 c------------------------------------------------------------------------------
2423       subroutine readrst
2424       implicit none
2425       include 'DIMENSIONS'
2426       include 'COMMON.CHAIN'
2427       include 'COMMON.IOUNITS'
2428       include 'COMMON.CONTROL'
2429       include 'COMMON.MD'
2430       include 'COMMON.QRESTR'
2431 #ifdef FIVEDIAG
2432        include 'COMMON.LAGRANGE.5diag'
2433 #else
2434        include 'COMMON.LAGRANGE'
2435 #endif
2436       integer i,j
2437       open(irest2,file=rest2name,status='unknown')
2438       read(irest2,*) totT,EK,potE,totE,t_bath
2439       totTafm=totT
2440       do i=0,2*nres-1
2441          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2442       enddo
2443       do i=0,2*nres-1
2444          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2445       enddo
2446       if(usampl) then
2447              read (irest2,*) iset
2448       endif
2449       close(irest2)
2450       return
2451       end
2452 c------------------------------------------------------------------------------
2453       subroutine read_fragments
2454       implicit none
2455       include 'DIMENSIONS'
2456 #ifdef MPI
2457       include 'mpif.h'
2458 #endif
2459       include 'COMMON.SETUP'
2460       include 'COMMON.CHAIN'
2461       include 'COMMON.IOUNITS'
2462       include 'COMMON.MD'
2463       include 'COMMON.QRESTR'
2464       include 'COMMON.CONTROL'
2465       integer i
2466       read(inp,*) nset,nfrag,npair,nfrag_back
2467       loc_qlike=(nfrag_back.lt.0)
2468       nfrag_back=iabs(nfrag_back)
2469       if(me.eq.king.or..not.out1file)
2470      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2471      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2472       do iset=1,nset
2473          read(inp,*) mset(iset)
2474        do i=1,nfrag
2475          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2476      &     qinfrag(i,iset)
2477          if(me.eq.king.or..not.out1file)
2478      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2479      &     ifrag(2,i,iset), qinfrag(i,iset)
2480        enddo
2481        do i=1,npair
2482         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2483      &    qinpair(i,iset)
2484         if(me.eq.king.or..not.out1file)
2485      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2486      &    ipair(2,i,iset), qinpair(i,iset)
2487        enddo
2488        if (loc_qlike) then
2489        do i=1,nfrag_back
2490         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2491      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2492      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2493      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2494         if(me.eq.king.or..not.out1file)
2495      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2496      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2497      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2498      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2499        enddo
2500        else
2501        do i=1,nfrag_back
2502         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2503      &     wfrag_back(3,i,iset),
2504      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2505         if(me.eq.king.or..not.out1file)
2506      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2507      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2508        enddo
2509        endif
2510       enddo
2511       return
2512       end
2513 C---------------------------------------------------------------------------
2514       subroutine read_afminp
2515       implicit none
2516       include 'DIMENSIONS'
2517 #ifdef MPI
2518       include 'mpif.h'
2519 #endif
2520       include 'COMMON.SETUP'
2521       include 'COMMON.CONTROL'
2522       include 'COMMON.CHAIN'
2523       include 'COMMON.IOUNITS'
2524       include 'COMMON.SBRIDGE'
2525       character*320 afmcard
2526       integer i
2527 c      print *, "wchodze"
2528       call card_concat(afmcard)
2529       call readi(afmcard,"BEG",afmbeg,0)
2530       call readi(afmcard,"END",afmend,0)
2531       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2532       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2533 c      print *,'FORCE=' ,forceAFMconst
2534 CCCC NOW PROPERTIES FOR AFM
2535        distafminit=0.0d0
2536        do i=1,3
2537         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2538        enddo
2539         distafminit=dsqrt(distafminit)
2540 c        print *,'initdist',distafminit
2541       return
2542       end
2543 c-------------------------------------------------------------------------------
2544       subroutine read_saxs_constr
2545       implicit none
2546       include 'DIMENSIONS'
2547 #ifdef MPI
2548       include 'mpif.h'
2549 #endif
2550       include 'COMMON.SETUP'
2551       include 'COMMON.CONTROL'
2552       include 'COMMON.SAXS'
2553       include 'COMMON.CHAIN'
2554       include 'COMMON.IOUNITS'
2555       include 'COMMON.SBRIDGE'
2556       double precision cm(3),cnorm
2557       integer i,j
2558 c      read(inp,*) nsaxs
2559       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2560       call flush(iout)
2561       if (saxs_mode.eq.0) then
2562 c SAXS distance distribution
2563       do i=1,nsaxs
2564         read(inp,*) distsaxs(i),Psaxs(i)
2565       enddo
2566       Cnorm = 0.0d0
2567       do i=1,nsaxs
2568         Cnorm = Cnorm + Psaxs(i)
2569       enddo
2570       write (iout,*) "Cnorm",Cnorm
2571       do i=1,nsaxs
2572         Psaxs(i)=Psaxs(i)/Cnorm
2573       enddo
2574       write (iout,*) "Normalized distance distribution from SAXS"
2575       do i=1,nsaxs
2576         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2577       enddo
2578       Wsaxs0=0.0d0
2579       do i=1,nsaxs
2580         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2581       enddo
2582       write (iout,*) "Wsaxs0",Wsaxs0
2583       else
2584 c SAXS "spheres".
2585       do i=1,nsaxs
2586         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2587       enddo
2588       do j=1,3
2589         cm(j)=0.0d0
2590       enddo
2591       do i=1,nsaxs
2592         do j=1,3
2593           cm(j)=cm(j)+Csaxs(j,i)
2594         enddo
2595       enddo
2596       do j=1,3
2597         cm(j)=cm(j)/nsaxs
2598       enddo
2599       do i=1,nsaxs
2600         do j=1,3
2601           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2602         enddo
2603       enddo
2604       write (iout,*) "SAXS sphere coordinates"
2605       do i=1,nsaxs
2606         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2607       enddo
2608       endif
2609       return
2610       end
2611
2612 c-------------------------------------------------------------------------------
2613       subroutine read_dist_constr
2614       implicit none
2615       include 'DIMENSIONS'
2616 #ifdef MPI
2617       include 'mpif.h'
2618 #endif
2619       include 'COMMON.SETUP'
2620       include 'COMMON.CONTROL'
2621       include 'COMMON.CHAIN'
2622       include 'COMMON.IOUNITS'
2623       include 'COMMON.SBRIDGE'
2624       include 'COMMON.INTERACT'
2625       integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2626       integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2627       double precision wfrag_(100),wpair_(1000)
2628       double precision ddjk,dist,dist_cut,fordepthmax
2629       character*5000 controlcard
2630       logical normalize,next
2631       integer restr_type
2632       double precision scal_bfac
2633       double precision xlink(4,0:4) /
2634 c           a          b       c     sigma
2635      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2636      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2637      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2638      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2639      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2640 c      print *, "WCHODZE"
2641       write (iout,*) "Calling read_dist_constr"
2642 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2643 c      call flush(iout)
2644       restr_on_coord=.false.
2645       next=.true.
2646
2647       npeak=0
2648       ipeak=0
2649       nhpb_peak=0
2650
2651       DO WHILE (next)
2652
2653       call card_concat(controlcard)
2654       next = index(controlcard,"NEXT").gt.0
2655       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2656       write (iout,*) "restr_type",restr_type
2657       call readi(controlcard,"NFRAG",nfrag_,0)
2658       call readi(controlcard,"NFRAG",nfrag_,0)
2659       call readi(controlcard,"NPAIR",npair_,0)
2660       call readi(controlcard,"NDIST",ndist_,0)
2661       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2662       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2663       if (restr_type.eq.10)
2664      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2665       if (restr_type.eq.12)
2666      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2667       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2668       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2669       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2670       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2671       normalize = index(controlcard,"NORMALIZE").gt.0
2672       write (iout,*) "WBOLTZD",wboltzd
2673       write (iout,*) "SCAL_PEAK",scal_peak
2674       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2675       write (iout,*) "IFRAG"
2676       do i=1,nfrag_
2677         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2678       enddo
2679       write (iout,*) "IPAIR"
2680       do i=1,npair_
2681         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2682       enddo
2683       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2684      & write (iout,*)
2685      &   "Distance restraints as generated from reference structure"
2686       do i=1,nfrag_
2687         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2688         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2689      &    ifrag_(2,i)=nstart_sup+nsup-1
2690 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2691 c        call flush(iout)
2692         if (wfrag_(i).eq.0.0d0) cycle
2693         do j=ifrag_(1,i),ifrag_(2,i)-1
2694           do k=j+1,ifrag_(2,i)
2695 c            write (iout,*) "j",j," k",k
2696             ddjk=dist(j,k)
2697             if (restr_type.eq.1) then
2698               nhpb=nhpb+1
2699               irestr_type(nhpb)=1
2700               ihpb(nhpb)=j
2701               jhpb(nhpb)=k
2702               dhpb(nhpb)=ddjk
2703               forcon(nhpb)=wfrag_(i)
2704             else if (constr_dist.eq.2) then
2705               if (ddjk.le.dist_cut) then
2706                 nhpb=nhpb+1
2707                 irestr_type(nhpb)=1
2708                 ihpb(nhpb)=j
2709                 jhpb(nhpb)=k
2710                 dhpb(nhpb)=ddjk
2711                 forcon(nhpb)=wfrag_(i)
2712               endif
2713             else if (restr_type.eq.3) then
2714               nhpb=nhpb+1
2715               irestr_type(nhpb)=1
2716               ihpb(nhpb)=j
2717               jhpb(nhpb)=k
2718               dhpb(nhpb)=ddjk
2719               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2720             endif
2721 #ifdef MPI
2722             if (.not.out1file .or. me.eq.king)
2723      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2724      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2725 #else
2726             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2727      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2728 #endif
2729           enddo
2730         enddo
2731       enddo
2732       do i=1,npair_
2733         if (wpair_(i).eq.0.0d0) cycle
2734         ii = ipair_(1,i)
2735         jj = ipair_(2,i)
2736         if (ii.gt.jj) then
2737           itemp=ii
2738           ii=jj
2739           jj=itemp
2740         endif
2741         do j=ifrag_(1,ii),ifrag_(2,ii)
2742           do k=ifrag_(1,jj),ifrag_(2,jj)
2743             ddjk=dist(j,k)
2744             if (restr_type.eq.1) then
2745               nhpb=nhpb+1
2746               irestr_type(nhpb)=1
2747               ihpb(nhpb)=j
2748               jhpb(nhpb)=k
2749               dhpb(nhpb)=ddjk
2750               forcon(nhpb)=wpair_(i)
2751             else if (constr_dist.eq.2) then
2752               if (ddjk.le.dist_cut) then
2753                 nhpb=nhpb+1
2754                 irestr_type(nhpb)=1
2755                 ihpb(nhpb)=j
2756                 jhpb(nhpb)=k
2757                 dhpb(nhpb)=ddjk
2758                 forcon(nhpb)=wpair_(i)
2759               endif
2760             else if (restr_type.eq.3) then
2761               nhpb=nhpb+1
2762               irestr_type(nhpb)=1
2763               ihpb(nhpb)=j
2764               jhpb(nhpb)=k
2765               dhpb(nhpb)=ddjk
2766               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2767             endif
2768 #ifdef MPI
2769             if (.not.out1file .or. me.eq.king)
2770      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2771      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2772 #else
2773             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2774      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2775 #endif
2776           enddo
2777         enddo
2778       enddo
2779
2780 c      print *,ndist_
2781       write (iout,*) "Distance restraints as read from input"
2782       do i=1,ndist_
2783         if (restr_type.eq.12) then
2784           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2785      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2786      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2787      &    fordepth_peak(nhpb_peak+1),npeak
2788 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2789 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2790 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2791 c     &    fordepth_peak(nhpb_peak+1),npeak
2792           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2793      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2794           nhpb_peak=nhpb_peak+1
2795           irestr_type_peak(nhpb_peak)=12
2796           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2797           ipeak(2,npeak)=i
2798 #ifdef MPI
2799           if (.not.out1file .or. me.eq.king)
2800      &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2801      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2802      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2803      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2804      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2805 #else
2806           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2807      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2808      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2809      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2810      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2811 #endif
2812           if (ibecarb_peak(nhpb_peak).eq.3) then
2813             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2814           else if (ibecarb_peak(nhpb_peak).eq.2) then
2815             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2816           else if (ibecarb_peak(nhpb_peak).eq.1) then
2817             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2818             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2819           endif
2820         else if (restr_type.eq.11) then
2821           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2822      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2823 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2824           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2825           nhpb=nhpb+1
2826           irestr_type(nhpb)=11
2827 #ifdef MPI
2828           if (.not.out1file .or. me.eq.king)
2829      &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2830      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2831      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2832 #else
2833           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2834      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2835      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2836 #endif
2837 c          if (ibecarb(nhpb).gt.0) then
2838 c            ihpb(nhpb)=ihpb(nhpb)+nres
2839 c            jhpb(nhpb)=jhpb(nhpb)+nres
2840 c          endif
2841          if (ibecarb(nhpb).eq.3) then
2842             ihpb(nhpb)=ihpb(nhpb)+nres
2843           else if (ibecarb(nhpb).eq.2) then
2844             ihpb(nhpb)=ihpb(nhpb)+nres
2845           else if (ibecarb(nhpb).eq.1) then
2846             ihpb(nhpb)=ihpb(nhpb)+nres
2847             jhpb(nhpb)=jhpb(nhpb)+nres
2848           endif
2849         else if (restr_type.eq.10) then
2850 c Cross-lonk Markov-like potential
2851           call card_concat(controlcard)
2852           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2853           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2854           ibecarb(nhpb+1)=0
2855           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2856           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2857           if (index(controlcard,"ZL").gt.0) then
2858             link_type=1
2859           else if (index(controlcard,"ADH").gt.0) then
2860             link_type=2
2861           else if (index(controlcard,"PDH").gt.0) then
2862             link_type=3
2863           else if (index(controlcard,"DSS").gt.0) then
2864             link_type=4
2865           else
2866             link_type=0
2867           endif
2868           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2869      &       xlink(1,link_type))
2870           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2871      &       xlink(2,link_type))
2872           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2873      &       xlink(3,link_type))
2874           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2875      &       xlink(4,link_type))
2876           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2877 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2878 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2879           if (forcon(nhpb+1).le.0.0d0 .or.
2880      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2881           nhpb=nhpb+1
2882           irestr_type(nhpb)=10
2883           if (ibecarb(nhpb).eq.3) then
2884             jhpb(nhpb)=jhpb(nhpb)+nres
2885           else if (ibecarb(nhpb).eq.2) then
2886             ihpb(nhpb)=ihpb(nhpb)+nres
2887           else if (ibecarb(nhpb).eq.1) then
2888             ihpb(nhpb)=ihpb(nhpb)+nres
2889             jhpb(nhpb)=jhpb(nhpb)+nres
2890           endif
2891 #ifdef MPI
2892           if (.not.out1file .or. me.eq.king)
2893      &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2894      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2895      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2896      &     irestr_type(nhpb)
2897 #else
2898           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2899      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2900      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2901      &     irestr_type(nhpb)
2902 #endif
2903         else
2904 C        print *,"in else"
2905           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2906      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2907           if (forcon(nhpb+1).gt.0.0d0) then
2908           nhpb=nhpb+1
2909           if (dhpb1(nhpb).eq.0.0d0) then
2910             irestr_type(nhpb)=1
2911           else
2912             irestr_type(nhpb)=2
2913           endif
2914           if (ibecarb(nhpb).eq.3) then
2915             jhpb(nhpb)=jhpb(nhpb)+nres
2916           else if (ibecarb(nhpb).eq.2) then
2917             ihpb(nhpb)=ihpb(nhpb)+nres
2918           else if (ibecarb(nhpb).eq.1) then
2919             ihpb(nhpb)=ihpb(nhpb)+nres
2920             jhpb(nhpb)=jhpb(nhpb)+nres
2921           endif
2922           if (dhpb(nhpb).eq.0.0d0)
2923      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2924         endif
2925 #ifdef MPI
2926           if (.not.out1file .or. me.eq.king)
2927      &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2928      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2929 #else
2930           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2931      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2932 #endif
2933         endif
2934 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2935 C        if (forcon(nhpb+1).gt.0.0d0) then
2936 C          nhpb=nhpb+1
2937 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2938       enddo
2939
2940       if (restr_type.eq.4) then
2941         write (iout,*) "The BFAC array"
2942         do i=nnt,nct
2943           write (iout,'(i5,f10.5)') i,bfac(i)
2944         enddo
2945         do i=nnt,nct
2946           if (itype(i).eq.ntyp1) cycle
2947           do j=nnt,i-1
2948             if (itype(j).eq.ntyp1) cycle
2949             if (itype(i).eq.10) then
2950               iiend=0
2951             else
2952               iiend=1
2953             endif
2954             if (itype(j).eq.10) then
2955               jjend=0
2956             else
2957               jjend=1
2958             endif
2959             kk=0
2960             do ii=0,iiend
2961             do jj=0,jjend
2962             nhpb=nhpb+1
2963             irestr_type(nhpb)=1
2964             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2965             irestr_type(nhpb)=1
2966             ibecarb(nhpb)=kk
2967             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2968             ihpb(nhpb)=i+nres*ii
2969             jhpb(nhpb)=j+nres*jj
2970             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2971 #ifdef MPI
2972             if (.not.out1file .or. me.eq.king) then
2973             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2974      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2975      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2976      &       irestr_type(nhpb)
2977             endif
2978 #else
2979             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2980      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2981      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2982      &       irestr_type(nhpb)
2983 #endif
2984             kk=kk+1
2985           enddo
2986           enddo
2987           enddo
2988         enddo
2989       endif
2990
2991       if (restr_type.eq.5) then
2992         restr_on_coord=.true.
2993         do i=nnt,nct
2994           if (itype(i).eq.ntyp1) cycle
2995           bfac(i)=(scal_bfac/bfac(i))**2
2996         enddo
2997       endif
2998
2999       ENDDO ! next
3000
3001       fordepthmax=0.0d0
3002       if (normalize) then
3003         do i=nss+1,nhpb
3004           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3005      &      fordepthmax=fordepth(i)
3006         enddo
3007         do i=nss+1,nhpb
3008           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3009         enddo
3010       endif
3011       return
3012       end
3013 c-------------------------------------------------------------------------------
3014       subroutine read_constr_homology
3015       implicit none
3016       include 'DIMENSIONS'
3017 #ifdef MPI
3018       include 'mpif.h'
3019 #endif
3020       include 'COMMON.SETUP'
3021       include 'COMMON.CONTROL'
3022       include 'COMMON.HOMOLOGY'
3023       include 'COMMON.CHAIN'
3024       include 'COMMON.IOUNITS'
3025       include 'COMMON.MD'
3026       include 'COMMON.QRESTR'
3027       include 'COMMON.GEO'
3028       include 'COMMON.INTERACT'
3029       include 'COMMON.NAMES'
3030 c
3031 c For new homol impl
3032 c
3033       include 'COMMON.VAR'
3034 c
3035
3036 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3037 c    &                 dist_cut
3038 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3039 c    &    sigma_odl_temp(maxres,maxres,max_template)
3040       character*2 kic2
3041       character*24 model_ki_dist, model_ki_angle
3042       character*500 controlcard
3043       integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
3044      & ik,iistart,nres_temp
3045       integer ilen
3046       external ilen
3047       logical liiflag,lfirst
3048       integer i01,i10
3049 c
3050 c     FP - Nov. 2014 Temporary specifications for new vars
3051 c
3052       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3053      &                 rescore3_tmp
3054       double precision, dimension (max_template,maxres) :: rescore
3055       double precision, dimension (max_template,maxres) :: rescore2
3056       double precision, dimension (max_template,maxres) :: rescore3
3057       double precision distal
3058       character*24 tpl_k_rescore
3059       character*256 pdbfile
3060 c -----------------------------------------------------------------
3061 c Reading multiple PDB ref structures and calculation of retraints
3062 c not using pre-computed ones stored in files model_ki_{dist,angle}
3063 c FP (Nov., 2014)
3064 c -----------------------------------------------------------------
3065 c
3066 c
3067 c Alternative: reading from input
3068       call card_concat(controlcard)
3069       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3070       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3071       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3072       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3073       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3074       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3075       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3076       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3077       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3078       if(.not.read2sigma.and.start_from_model) then
3079           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3080      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3081           start_from_model=.false.
3082           iranconf=(indpdb.le.0)
3083       endif
3084       if(start_from_model .and. (me.eq.king .or. .not. out1file))
3085      &    write(iout,*) 'START_FROM_MODELS is ON'
3086 c      if(start_from_model .and. rest) then
3087 c        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3088 c         write(iout,*) 'START_FROM_MODELS is OFF'
3089 c         write(iout,*) 'remove restart keyword from input'
3090 c        endif
3091 c      endif
3092       if (start_from_model) nmodel_start=constr_homology
3093       if (homol_nset.gt.1)then
3094          call card_concat(controlcard)
3095          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3096          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3097           write(iout,*) "iset homology_weight "
3098           do i=1,homol_nset
3099            write(iout,*) i,waga_homology(i)
3100           enddo
3101          endif
3102          iset=mod(kolor,homol_nset)+1
3103       else
3104        iset=1
3105        waga_homology(1)=1.0
3106       endif
3107
3108 cd      write (iout,*) "nnt",nnt," nct",nct
3109 cd      call flush(iout)
3110
3111
3112       lim_odl=0
3113       lim_dih=0
3114 c
3115 c      write(iout,*) 'nnt=',nnt,'nct=',nct
3116 c
3117       do i = nnt,nct
3118         do k=1,constr_homology
3119          idomain(k,i)=0
3120         enddo
3121       enddo
3122
3123       ii=0
3124       do i = nnt,nct-2
3125         do j=i+2,nct
3126         ii=ii+1
3127         ii_in_use(ii)=0
3128         enddo
3129       enddo
3130
3131       if (read_homol_frag) then
3132         call read_klapaucjusz
3133       else
3134
3135       do k=1,constr_homology
3136
3137         read(inp,'(a)') pdbfile
3138         pdbfiles_chomo(k)=pdbfile
3139         if(me.eq.king .or. .not. out1file)
3140      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3141      &  pdbfile(:ilen(pdbfile))
3142         open(ipdbin,file=pdbfile,status='old',err=33)
3143         goto 34
3144   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3145      &  pdbfile(:ilen(pdbfile))
3146         stop
3147   34    continue
3148 c        print *,'Begin reading pdb data'
3149 c
3150 c Files containing res sim or local scores (former containing sigmas)
3151 c
3152
3153         write(kic2,'(bz,i2.2)') k
3154
3155         tpl_k_rescore="template"//kic2//".sco"
3156
3157         unres_pdb=.false.
3158         nres_temp=nres
3159         if (read2sigma) then
3160           call readpdb_template(k)
3161         else
3162           call readpdb
3163         endif
3164         nres_chomo(k)=nres
3165         nres=nres_temp
3166 c
3167 c     Distance restraints
3168 c
3169 c          ... --> odl(k,ii)
3170 C Copy the coordinates from reference coordinates (?)
3171         do i=1,2*nres_chomo(k)
3172           do j=1,3
3173             c(j,i)=cref(j,i)
3174 c           write (iout,*) "c(",j,i,") =",c(j,i)
3175           enddo
3176         enddo
3177 c
3178 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3179 c
3180 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3181           open (ientin,file=tpl_k_rescore,status='old')
3182           if (nnt.gt.1) rescore(k,1)=0.0d0
3183           do irec=nnt,nct ! loop for reading res sim
3184             if (read2sigma) then
3185              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3186      &                                rescore3_tmp,idomain_tmp
3187              i_tmp=i_tmp+nnt-1
3188              idomain(k,i_tmp)=idomain_tmp
3189              rescore(k,i_tmp)=rescore_tmp
3190              rescore2(k,i_tmp)=rescore2_tmp
3191              rescore3(k,i_tmp)=rescore3_tmp
3192              if (.not. out1file .or. me.eq.king)
3193      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3194      &                      i_tmp,rescore2_tmp,rescore_tmp,
3195      &                                rescore3_tmp,idomain_tmp
3196             else
3197              idomain(k,irec)=1
3198              read (ientin,*,end=1401) rescore_tmp
3199
3200 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3201              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3202 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3203             endif
3204           enddo
3205  1401   continue
3206         close (ientin)
3207         if (waga_dist.ne.0.0d0) then
3208           ii=0
3209           do i = nnt,nct-2
3210             do j=i+2,nct
3211
3212               x12=c(1,i)-c(1,j)
3213               y12=c(2,i)-c(2,j)
3214               z12=c(3,i)-c(3,j)
3215               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3216 c              write (iout,*) k,i,j,distal,dist2_cut
3217
3218             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3219      &            .and. distal.le.dist2_cut ) then
3220
3221               ii=ii+1
3222               ii_in_use(ii)=1
3223               l_homo(k,ii)=.true.
3224
3225 c             write (iout,*) "k",k
3226 c             write (iout,*) "i",i," j",j," constr_homology",
3227 c    &                       constr_homology
3228               ires_homo(ii)=i
3229               jres_homo(ii)=j
3230               odl(k,ii)=distal
3231               if (read2sigma) then
3232                 sigma_odl(k,ii)=0
3233                 do ik=i,j
3234                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3235                 enddo
3236                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3237                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3238      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3239               else
3240                 if (odl(k,ii).le.dist_cut) then
3241                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3242                 else
3243 #ifdef OLDSIGMA
3244                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3245      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3246 #else
3247                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3248      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3249 #endif
3250                 endif
3251               endif
3252               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3253             else
3254               ii=ii+1
3255               l_homo(k,ii)=.false.
3256             endif
3257             enddo
3258           enddo
3259         lim_odl=ii
3260         endif
3261 c        write (iout,*) "Distance restraints set"
3262 c        call flush(iout)
3263 c
3264 c     Theta, dihedral and SC retraints
3265 c
3266         if (waga_angle.gt.0.0d0) then
3267 c         open (ientin,file=tpl_k_sigma_dih,status='old')
3268 c         do irec=1,maxres-3 ! loop for reading sigma_dih
3269 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3270 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3271 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3272 c    &                            sigma_dih(k,i+nnt-1)
3273 c         enddo
3274 c1402   continue
3275 c         close (ientin)
3276           do i = nnt+3,nct
3277             if (idomain(k,i).eq.0) then
3278                sigma_dih(k,i)=0.0
3279                cycle
3280             endif
3281             dih(k,i)=phiref(i) ! right?
3282 c           read (ientin,*) sigma_dih(k,i) ! original variant
3283 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
3284 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3285 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3286 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
3287 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
3288
3289             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3290      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3291 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3292 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
3293 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3294 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3295 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
3296 c   Instead of res sim other local measure of b/b str reliability possible
3297             if (sigma_dih(k,i).ne.0)
3298      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3299 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3300           enddo
3301           lim_dih=nct-nnt-2
3302         endif
3303 c        write (iout,*) "Dihedral angle restraints set"
3304 c        call flush(iout)
3305
3306         if (waga_theta.gt.0.0d0) then
3307 c         open (ientin,file=tpl_k_sigma_theta,status='old')
3308 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3309 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3310 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3311 c    &                              sigma_theta(k,i+nnt-1)
3312 c         enddo
3313 c1403   continue
3314 c         close (ientin)
3315
3316           do i = nnt+2,nct ! right? without parallel.
3317 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
3318 c         do i=ithet_start,ithet_end ! with FG parallel.
3319              if (idomain(k,i).eq.0) then
3320               sigma_theta(k,i)=0.0
3321               cycle
3322              endif
3323              thetatpl(k,i)=thetaref(i)
3324 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3325 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
3326 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3327 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
3328 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
3329              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3330      &                        rescore(k,i-2))/3.0
3331 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3332              if (sigma_theta(k,i).ne.0)
3333      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3334
3335 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3336 c                             rescore(k,i-2) !  right expression ?
3337 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3338           enddo
3339         endif
3340 c        write (iout,*) "Angle restraints set"
3341 c        call flush(iout)
3342
3343         if (waga_d.gt.0.0d0) then
3344 c       open (ientin,file=tpl_k_sigma_d,status='old')
3345 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3346 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3347 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3348 c    &                          sigma_d(k,i+nnt-1)
3349 c         enddo
3350 c1404   continue
3351
3352           do i = nnt,nct ! right? without parallel.
3353 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3354 c         do i=loc_start,loc_end ! with FG parallel.
3355                if (itype(i).eq.10) cycle
3356                if (idomain(k,i).eq.0 ) then
3357                   sigma_d(k,i)=0.0
3358                   cycle
3359                endif
3360                xxtpl(k,i)=xxref(i)
3361                yytpl(k,i)=yyref(i)
3362                zztpl(k,i)=zzref(i)
3363 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3364 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3365 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3366 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3367                sigma_d(k,i)=rescore3(k,i) !  right expression ?
3368                if (sigma_d(k,i).ne.0)
3369      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3370
3371 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3372 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3373 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3374           enddo
3375         endif
3376       enddo
3377 c      write (iout,*) "SC restraints set"
3378 c      call flush(iout)
3379 c
3380 c remove distance restraints not used in any model from the list
3381 c shift data in all arrays
3382 c
3383 c      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3384       if (waga_dist.ne.0.0d0) then
3385         ii=0
3386         liiflag=.true.
3387         lfirst=.true.
3388         do i=nnt,nct-2
3389          do j=i+2,nct
3390           ii=ii+1
3391 c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3392 c     &            .and. distal.le.dist2_cut ) then
3393 c          write (iout,*) "i",i," j",j," ii",ii
3394 c          call flush(iout)
3395           if (ii_in_use(ii).eq.0.and.liiflag.or.
3396      &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3397             liiflag=.false.
3398             i10=ii
3399             if (lfirst) then
3400               lfirst=.false.
3401               iistart=ii
3402             else
3403               if(i10.eq.lim_odl) i10=i10+1
3404               do ki=0,i10-i01-1
3405                ires_homo(iistart+ki)=ires_homo(ki+i01)
3406                jres_homo(iistart+ki)=jres_homo(ki+i01)
3407                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3408                do k=1,constr_homology
3409                 odl(k,iistart+ki)=odl(k,ki+i01)
3410                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3411                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3412                enddo
3413               enddo
3414               iistart=iistart+i10-i01
3415             endif
3416           endif
3417           if (ii_in_use(ii).ne.0.and..not.liiflag) then
3418              i01=ii
3419              liiflag=.true.
3420           endif
3421          enddo
3422         enddo
3423         lim_odl=iistart-1
3424       endif
3425 c      write (iout,*) "Removing distances completed"
3426 c      call flush(iout)
3427       endif ! .not. klapaucjusz
3428
3429       if (constr_homology.gt.0) call homology_partition
3430 c      write (iout,*) "After homology_partition"
3431 c      call flush(iout)
3432       if (constr_homology.gt.0) call init_int_table
3433 c      write (iout,*) "After init_int_table"
3434 c      call flush(iout)
3435 c      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3436 c      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3437 c
3438 c Print restraints
3439 c
3440       if (.not.out_template_restr) return
3441 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3442       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3443        write (iout,*) "Distance restraints from templates"
3444        do ii=1,lim_odl
3445        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3446      &  ii,ires_homo(ii),jres_homo(ii),
3447      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3448      &  ki=1,constr_homology)
3449        enddo
3450        write (iout,*) "Dihedral angle restraints from templates"
3451        do i=nnt+3,nct
3452         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3453      &      (rad2deg*dih(ki,i),
3454      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3455        enddo
3456        write (iout,*) "Virtual-bond angle restraints from templates"
3457        do i=nnt+2,nct
3458         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3459      &      (rad2deg*thetatpl(ki,i),
3460      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3461        enddo
3462        write (iout,*) "SC restraints from templates"
3463        do i=nnt,nct
3464         write(iout,'(i5,100(4f8.2,4x))') i,
3465      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3466      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3467        enddo
3468       endif
3469 c -----------------------------------------------------------------
3470       return
3471       end
3472 c----------------------------------------------------------------------
3473 #ifdef WINIFL
3474       subroutine flush(iu)
3475       return
3476       end
3477 #endif
3478 #ifdef AIX
3479       subroutine flush(iu)
3480       call flush_(iu)
3481       return
3482       end
3483 #endif
3484 c------------------------------------------------------------------------------
3485       subroutine copy_to_tmp(source)
3486       implicit none
3487       include "DIMENSIONS"
3488       include "COMMON.IOUNITS"
3489       character*(*) source
3490       character* 256 tmpfile
3491       integer ilen
3492       external ilen
3493       logical ex
3494       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3495       inquire(file=tmpfile,exist=ex)
3496       if (ex) then
3497         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3498      &   " to temporary directory..."
3499         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3500         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3501       endif
3502       return
3503       end
3504 c------------------------------------------------------------------------------
3505       subroutine move_from_tmp(source)
3506       implicit none
3507       include "DIMENSIONS"
3508       include "COMMON.IOUNITS"
3509       character*(*) source
3510       integer ilen
3511       external ilen
3512       write (*,*) "Moving ",source(:ilen(source)),
3513      & " from temporary directory to working directory"
3514       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3515       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3516       return
3517       end
3518 c------------------------------------------------------------------------------
3519       subroutine random_init(seed)
3520 C
3521 C Initialize random number generator
3522 C
3523       implicit none
3524       include 'DIMENSIONS'
3525 #ifdef MPI
3526       include 'mpif.h'
3527       logical OKRandom, prng_restart
3528       real*8  r1
3529       integer iseed_array(4)
3530       integer error_msg,ierr
3531 #endif
3532       include 'COMMON.IOUNITS'
3533       include 'COMMON.TIME1'
3534       include 'COMMON.THREAD'
3535       include 'COMMON.SBRIDGE'
3536       include 'COMMON.CONTROL'
3537       include 'COMMON.MCM'
3538       include 'COMMON.MAP'
3539       include 'COMMON.HEADER'
3540       include 'COMMON.CSA'
3541       include 'COMMON.CHAIN'
3542       include 'COMMON.MUCA'
3543       include 'COMMON.MD'
3544       include 'COMMON.FFIELD'
3545       include 'COMMON.SETUP'
3546       integer i,iseed
3547       double precision seed,ran_number
3548       iseed=-dint(dabs(seed))
3549       if (iseed.eq.0) then
3550         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3551      &    'Random seed undefined. The program will stop.'
3552         write (*,'(/80(1h*)/20x,a/80(1h*))')
3553      &    'Random seed undefined. The program will stop.'
3554 #ifdef MPI
3555         call mpi_finalize(mpi_comm_world,ierr)
3556 #endif
3557         stop 'Bad random seed.'
3558       endif
3559 #ifdef MPI
3560       if (fg_rank.eq.0) then
3561       seed=seed*(me+1)+1
3562 #ifdef AMD64
3563       if(me.eq.king)
3564      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3565       OKRandom = prng_restart(me,iseed)
3566 #else
3567       do i=1,4
3568        tmp=65536.0d0**(4-i)
3569        iseed_array(i) = dint(seed/tmp)
3570        seed=seed-iseed_array(i)*tmp
3571       enddo
3572       if(me.eq.king)
3573      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3574      &                 (iseed_array(i),i=1,4)
3575       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3576      &                 (iseed_array(i),i=1,4)
3577       OKRandom = prng_restart(me,iseed_array)
3578 #endif
3579       if (OKRandom) then
3580 c        r1 = prng_next(me)
3581         r1=ran_number(0.0D0,1.0D0)
3582         if(me.eq.king)
3583      &   write (iout,*) 'ran_num',r1
3584         if (r1.lt.0.0d0) OKRandom=.false.
3585       endif
3586       if (.not.OKRandom) then
3587         write (iout,*) 'PRNG IS NOT WORKING!!!'
3588         print *,'PRNG IS NOT WORKING!!!'
3589         if (me.eq.0) then
3590          call flush(iout)
3591          call mpi_abort(mpi_comm_world,error_msg,ierr)
3592          stop
3593         else
3594          write (iout,*) 'too many processors for parallel prng'
3595          write (*,*) 'too many processors for parallel prng'
3596          call flush(iout)
3597          stop
3598         endif
3599       endif
3600       endif
3601 #else
3602       call vrndst(iseed)
3603       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3604 #endif
3605       return
3606       end
3607 c----------------------------------------------------------------------
3608       subroutine read_klapaucjusz
3609       implicit none
3610       include 'DIMENSIONS'
3611 #ifdef MPI
3612       include 'mpif.h'
3613 #endif
3614       include 'COMMON.SETUP'
3615       include 'COMMON.CONTROL'
3616       include 'COMMON.HOMOLOGY'
3617       include 'COMMON.CHAIN'
3618       include 'COMMON.IOUNITS'
3619       include 'COMMON.MD'
3620       include 'COMMON.GEO'
3621       include 'COMMON.INTERACT'
3622       include 'COMMON.NAMES'
3623       character*256 fragfile
3624       integer ninclust(maxclust),inclust(max_template,maxclust),
3625      &  nresclust(maxclust),iresclust(maxres,maxclust),nclust
3626
3627       character*2 kic2
3628       character*24 model_ki_dist, model_ki_angle
3629       character*500 controlcard
3630       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3631      & ik,ll,ii,kk,iistart,iishift,lim_xx
3632       double precision distal
3633       logical lprn /.true./
3634       integer nres_temp
3635       integer ilen
3636       external ilen
3637       logical liiflag
3638 c
3639 c
3640       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3641       double precision, dimension (max_template,maxres) :: rescore
3642       double precision, dimension (max_template,maxres) :: rescore2
3643       character*24 tpl_k_rescore
3644       character*256 pdbfile
3645
3646 c
3647 c For new homol impl
3648 c
3649       include 'COMMON.VAR'
3650 c
3651       call getenv("FRAGFILE",fragfile)
3652       open(ientin,file=fragfile,status="old",err=10)
3653       read(ientin,*) constr_homology,nclust
3654       l_homo = .false.
3655       sigma_theta=0.0
3656       sigma_d=0.0
3657       sigma_dih=0.0
3658 c Read pdb files
3659       do k=1,constr_homology
3660         read(ientin,'(a)') pdbfile
3661         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3662      &  pdbfile(:ilen(pdbfile))
3663         pdbfiles_chomo(k)=pdbfile
3664         open(ipdbin,file=pdbfile,status='old',err=33)
3665         goto 34
3666   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3667      &  pdbfile(:ilen(pdbfile))
3668         stop
3669   34    continue
3670         unres_pdb=.false.
3671         nres_temp=nres
3672         call readpdb_template(k)
3673         nres_chomo(k)=nres
3674         nres=nres_temp
3675         do i=1,nres
3676           rescore(k,i)=0.2d0
3677           rescore2(k,i)=1.0d0
3678         enddo
3679       enddo
3680 c Read clusters
3681       do i=1,nclust
3682         read(ientin,*) ninclust(i),nresclust(i)
3683         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3684         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3685       enddo
3686 c
3687 c Loop over clusters
3688 c
3689       do l=1,nclust
3690         do ll = 1,ninclust(l)
3691
3692         k = inclust(ll,l)
3693         do i=1,nres
3694           idomain(k,i)=0
3695         enddo
3696         do i=1,nresclust(l)
3697           if (nnt.gt.1)  then
3698             idomain(k,iresclust(i,l)+1) = 1
3699           else
3700             idomain(k,iresclust(i,l)) = 1
3701           endif
3702         enddo
3703 c
3704 c     Distance restraints
3705 c
3706 c          ... --> odl(k,ii)
3707 C Copy the coordinates from reference coordinates (?)
3708         nres_temp=nres
3709         nres=nres_chomo(k)
3710         do i=1,2*nres
3711           do j=1,3
3712             c(j,i)=chomo(j,i,k)
3713 c           write (iout,*) "c(",j,i,") =",c(j,i)
3714           enddo
3715         enddo
3716         call int_from_cart(.true.,.false.)
3717         call sc_loc_geom(.false.)
3718         do i=1,nres
3719           thetaref(i)=theta(i)
3720           phiref(i)=phi(i)
3721         enddo
3722         nres=nres_temp
3723         if (waga_dist.ne.0.0d0) then
3724           ii=0
3725           do i = nnt,nct-2
3726             do j=i+2,nct
3727
3728               x12=c(1,i)-c(1,j)
3729               y12=c(2,i)-c(2,j)
3730               z12=c(3,i)-c(3,j)
3731               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3732 c              write (iout,*) k,i,j,distal,dist2_cut
3733
3734             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3735      &            .and. distal.le.dist2_cut ) then
3736
3737               ii=ii+1
3738               ii_in_use(ii)=1
3739               l_homo(k,ii)=.true.
3740
3741 c             write (iout,*) "k",k
3742 c             write (iout,*) "i",i," j",j," constr_homology",
3743 c    &                       constr_homology
3744               ires_homo(ii)=i
3745               jres_homo(ii)=j
3746               odl(k,ii)=distal
3747               if (read2sigma) then
3748                 sigma_odl(k,ii)=0
3749                 do ik=i,j
3750                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3751                 enddo
3752                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3753                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3754      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3755               else
3756                 if (odl(k,ii).le.dist_cut) then
3757                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3758                 else
3759 #ifdef OLDSIGMA
3760                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3761      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3762 #else
3763                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3764      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3765 #endif
3766                 endif
3767               endif
3768               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3769             else
3770               ii=ii+1
3771 c              l_homo(k,ii)=.false.
3772             endif
3773             enddo
3774           enddo
3775         lim_odl=ii
3776         endif
3777 c
3778 c     Theta, dihedral and SC retraints
3779 c
3780         if (waga_angle.gt.0.0d0) then
3781           do i = nnt+3,nct
3782             if (idomain(k,i).eq.0) then
3783 c               sigma_dih(k,i)=0.0
3784                cycle
3785             endif
3786             dih(k,i)=phiref(i)
3787             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3788      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3789 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3790 c     &       " sigma_dihed",sigma_dih(k,i)
3791             if (sigma_dih(k,i).ne.0)
3792      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3793           enddo
3794           lim_dih=nct-nnt-2
3795         endif
3796
3797         if (waga_theta.gt.0.0d0) then
3798           do i = nnt+2,nct
3799              if (idomain(k,i).eq.0) then
3800 c              sigma_theta(k,i)=0.0
3801               cycle
3802              endif
3803              thetatpl(k,i)=thetaref(i)
3804              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3805      &                        rescore(k,i-2))/3.0
3806              if (sigma_theta(k,i).ne.0)
3807      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3808           enddo
3809         endif
3810
3811         if (waga_d.gt.0.0d0) then
3812           do i = nnt,nct
3813                if (itype(i).eq.10) cycle
3814                if (idomain(k,i).eq.0 ) then
3815 c                  sigma_d(k,i)=0.0
3816                   cycle
3817                endif
3818                xxtpl(k,i)=xxref(i)
3819                yytpl(k,i)=yyref(i)
3820                zztpl(k,i)=zzref(i)
3821                sigma_d(k,i)=rescore(k,i)
3822                if (sigma_d(k,i).ne.0)
3823      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3824                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3825           enddo
3826         endif
3827       enddo ! l
3828       enddo ! ll
3829 c
3830 c remove distance restraints not used in any model from the list
3831 c shift data in all arrays
3832 c
3833       if (waga_dist.ne.0.0d0) then
3834         ii=0
3835         liiflag=.true.
3836         do i=nnt,nct-2
3837          do j=i+2,nct
3838           ii=ii+1
3839           if (ii_in_use(ii).eq.0.and.liiflag) then
3840             liiflag=.false.
3841             iistart=ii
3842           endif
3843           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3844      &                   .not.liiflag.and.ii.eq.lim_odl) then
3845              if (ii.eq.lim_odl) then
3846               iishift=ii-iistart+1
3847              else
3848               iishift=ii-iistart
3849              endif
3850              liiflag=.true.
3851              do ki=iistart,lim_odl-iishift
3852               ires_homo(ki)=ires_homo(ki+iishift)
3853               jres_homo(ki)=jres_homo(ki+iishift)
3854               ii_in_use(ki)=ii_in_use(ki+iishift)
3855               do k=1,constr_homology
3856                odl(k,ki)=odl(k,ki+iishift)
3857                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3858                l_homo(k,ki)=l_homo(k,ki+iishift)
3859               enddo
3860              enddo
3861              ii=ii-iishift
3862              lim_odl=lim_odl-iishift
3863           endif
3864          enddo
3865         enddo
3866       endif
3867
3868       return
3869    10 stop "Error in fragment file"
3870       end