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