2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
16 C Read control parameters for energy minimzation if required
17 if (minim) call read_minim
18 C Read MCM control parameters if required
19 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31 csa if (modecalc.eq.8) then
32 csa inquire (file="fort.40",exist=file_exist)
33 csa if (.not.file_exist) call csaread
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
48 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
76 csa include 'COMMON.CSA'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
90 read (INP,'(a)') titel
91 call card_concat(controlcard)
92 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94 call reada(controlcard,'SEED',seed,0.0D0)
95 call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97 read_cart=index(controlcard,'READ_CART').gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
100 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107 call reada(controlcard,'DRMS',drms,0.1D0)
108 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
110 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
111 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
112 write (iout,'(a,f10.1)')'DRMS = ',drms
113 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
114 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
116 call readi(controlcard,'NZ_START',nz_start,0)
117 call readi(controlcard,'NZ_END',nz_end,0)
118 call readi(controlcard,'IZ_SC',iz_sc,0)
120 safety = 60.0d0*safety
123 call reada(controlcard,"T_BATH",t_bath,300.0d0)
124 minim=(index(controlcard,'MINIMIZE').gt.0)
125 dccart=(index(controlcard,'CART').gt.0)
126 overlapsc=(index(controlcard,'OVERLAP').gt.0)
127 overlapsc=.not.overlapsc
128 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129 searchsc=.not.searchsc
130 sideadd=(index(controlcard,'SIDEADD').gt.0)
131 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132 outpdb=(index(controlcard,'PDBOUT').gt.0)
133 outmol2=(index(controlcard,'MOL2OUT').gt.0)
134 pdbref=(index(controlcard,'PDBREF').gt.0)
135 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136 indpdb=index(controlcard,'PDBSTART')
137 extconf=(index(controlcard,'EXTCONF').gt.0)
138 call readi(controlcard,'IPRINT',iprint,0)
139 call readi(controlcard,'MAXGEN',maxgen,10000)
140 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141 call readi(controlcard,"KDIAG",kdiag,0)
142 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144 & write (iout,*) "RESCALE_MODE",rescale_mode
145 split_ene=index(controlcard,'SPLIT_ENE').gt.0
146 if (index(controlcard,'REGULAR').gt.0.0D0) then
147 call reada(controlcard,'WEIDIS',weidis,0.1D0)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
180 write(*,*) "CSA not supported in this version"
183 crc else if (index(controlcard,'ZSCORE').gt.0) then
185 crc ZSCORE is rm from UNRES, modecalc=9 is available
188 cfcm else if (index(controlcard,'MCMF').gt.0) then
190 else if (index(controlcard,'SOFTREG').gt.0) then
192 else if (index(controlcard,'CHECK_BOND').gt.0) then
194 else if (index(controlcard,'TEST').gt.0) then
196 else if (index(controlcard,'MD').gt.0) then
198 else if (index(controlcard,'RE ').gt.0) then
202 lmuca=index(controlcard,'MUCA').gt.0
203 call readi(controlcard,'MUCADYN',mucadyn,0)
204 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
205 if (lmuca .and. (me.eq.king .or. .not.out1file ))
207 write (iout,*) 'MUCADYN=',mucadyn
208 write (iout,*) 'MUCASMOOTH=',muca_smooth
211 iscode=index(controlcard,'ONE_LETTER')
212 indphi=index(controlcard,'PHI')
213 indback=index(controlcard,'BACK')
214 iranconf=index(controlcard,'RAND_CONF')
215 i2ndstr=index(controlcard,'USE_SEC_PRED')
216 gradout=index(controlcard,'GRADOUT').gt.0
217 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
219 if(me.eq.king.or..not.out1file)
220 & write (iout,'(2a)') diagmeth(kdiag),
221 & ' routine used to diagonalize matrices.'
224 c--------------------------------------------------------------------------
225 subroutine read_REMDpar
229 implicit real*8 (a-h,o-z)
231 include 'COMMON.IOUNITS'
232 include 'COMMON.TIME1'
235 include 'COMMON.LANGEVIN'
237 include 'COMMON.LANGEVIN.lang0'
239 include 'COMMON.INTERACT'
240 include 'COMMON.NAMES'
242 include 'COMMON.REMD'
243 include 'COMMON.CONTROL'
244 include 'COMMON.SETUP'
246 character*320 controlcard
247 character*3200 controlcard1
248 integer iremd_m_total
250 if(me.eq.king.or..not.out1file)
251 & write (iout,*) "REMD setup"
253 call card_concat(controlcard)
254 call readi(controlcard,"NREP",nrep,3)
255 call readi(controlcard,"NSTEX",nstex,1000)
256 call reada(controlcard,"RETMIN",retmin,10.0d0)
257 call reada(controlcard,"RETMAX",retmax,1000.0d0)
258 mremdsync=(index(controlcard,'SYNC').gt.0)
259 call readi(controlcard,"NSYN",i_sync_step,100)
260 restart1file=(index(controlcard,'REST1FILE').gt.0)
261 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
262 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
263 if(max_cache_traj_use.gt.max_cache_traj)
264 & max_cache_traj_use=max_cache_traj
265 if(me.eq.king.or..not.out1file) then
266 cd if (traj1file) then
267 crc caching is in testing - NTWX is not ignored
268 cd write (iout,*) "NTWX value is ignored"
269 cd write (iout,*) " trajectory is stored to one file by master"
270 cd write (iout,*) " before exchange at NSTEX intervals"
272 write (iout,*) "NREP= ",nrep
273 write (iout,*) "NSTEX= ",nstex
274 write (iout,*) "SYNC= ",mremdsync
275 write (iout,*) "NSYN= ",i_sync_step
276 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
279 t_exchange_only=(index(controlcard,'TONLY').gt.0)
280 call readi(controlcard,"HREMD",hremd,0)
281 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
282 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
284 if(usampl.and.hremd.gt.0) then
286 & "========== ERROR: USAMPL and HREMD cannot be used together"
288 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
295 if (index(controlcard,'TLIST').gt.0) then
297 call card_concat(controlcard1)
298 read(controlcard1,*) (remd_t(i),i=1,nrep)
299 if(me.eq.king.or..not.out1file)
300 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
303 if (index(controlcard,'MLIST').gt.0) then
305 call card_concat(controlcard1)
306 read(controlcard1,*) (remd_m(i),i=1,nrep)
307 if(me.eq.king.or..not.out1file) then
308 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
311 iremd_m_total=iremd_m_total+remd_m(i)
314 write (iout,*) 'Total number of replicas ',
315 & iremd_m_total*hremd
317 write (iout,*) 'Total number of replicas ',iremd_m_total
321 if(me.eq.king.or..not.out1file)
322 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
325 c--------------------------------------------------------------------------
326 subroutine read_MDpar
330 implicit real*8 (a-h,o-z)
332 include 'COMMON.IOUNITS'
333 include 'COMMON.TIME1'
336 include 'COMMON.LANGEVIN'
338 include 'COMMON.LANGEVIN.lang0'
340 include 'COMMON.INTERACT'
341 include 'COMMON.NAMES'
343 include 'COMMON.SETUP'
344 include 'COMMON.CONTROL'
345 include 'COMMON.SPLITELE'
347 character*320 controlcard
349 call card_concat(controlcard)
350 call readi(controlcard,"NSTEP",n_timestep,1000000)
351 call readi(controlcard,"NTWE",ntwe,100)
352 call readi(controlcard,"NTWX",ntwx,1000)
353 call reada(controlcard,"DT",d_time,1.0d-1)
354 call reada(controlcard,"DVMAX",dvmax,2.0d1)
355 call reada(controlcard,"DAMAX",damax,1.0d1)
356 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
357 call readi(controlcard,"LANG",lang,0)
358 RESPA = index(controlcard,"RESPA") .gt. 0
359 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
360 ntime_split0=ntime_split
361 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
362 ntime_split0=ntime_split
363 call reada(controlcard,"R_CUT",r_cut,2.0d0)
364 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
365 rest = index(controlcard,"REST").gt.0
366 tbf = index(controlcard,"TBF").gt.0
367 call readi(controlcard,"HMC",hmc,0)
368 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
369 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
370 tnh = index(controlcard,"NOSEHOOVER96").gt.0
371 if (RESPA.and.tnh)then
372 xiresp = index(controlcard,"XIRESP").gt.0
374 call reada(controlcard,"Q_NP",Q_np,0.1d0)
375 usampl = index(controlcard,"USAMPL").gt.0
376 mdpdb = index(controlcard,"MDPDB").gt.0
377 call reada(controlcard,"T_BATH",t_bath,300.0d0)
378 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
379 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
380 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
381 if (count_reset_moment.eq.0) count_reset_moment=1000000000
382 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
383 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
384 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
385 if (count_reset_vel.eq.0) count_reset_vel=1000000000
386 large = index(controlcard,"LARGE").gt.0
387 print_compon = index(controlcard,"PRINT_COMPON").gt.0
388 rattle = index(controlcard,"RATTLE").gt.0
389 c if performing umbrella sampling, fragments constrained are read from the fragment file
395 if(me.eq.king.or..not.out1file) then
397 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
399 write (iout,'(a)') "The units are:"
400 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
401 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
402 & " acceleration: angstrom/(48.9 fs)**2"
403 write (iout,'(a)') "energy: kcal/mol, temperature: K"
405 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
406 write (iout,'(a60,f10.5,a)')
407 & "Initial time step of numerical integration:",d_time,
409 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
411 write (iout,'(2a,i4,a)')
412 & "A-MTS algorithm used; initial time step for fast-varying",
413 & " short-range forces split into",ntime_split," steps."
414 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
415 & r_cut," lambda",rlamb
417 write (iout,'(2a,f10.5)')
418 & "Maximum acceleration threshold to reduce the time step",
419 & "/increase split number:",damax
420 write (iout,'(2a,f10.5)')
421 & "Maximum predicted energy drift to reduce the timestep",
422 & "/increase split number:",edriftmax
423 write (iout,'(a60,f10.5)')
424 & "Maximum velocity threshold to reduce velocities:",dvmax
425 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
426 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
427 if (rattle) write (iout,'(a60)')
428 & "Rattle algorithm used to constrain the virtual bonds"
432 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
433 call reada(controlcard,"RWAT",rwat,1.4d0)
434 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
435 surfarea=index(controlcard,"SURFAREA").gt.0
436 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
437 if(me.eq.king.or..not.out1file)then
438 write (iout,'(/a,$)') "Langevin dynamics calculation"
441 & " with direct integration of Langevin equations"
442 else if (lang.eq.2) then
443 write (iout,'(a/)') " with TINKER stochasic MD integrator"
444 else if (lang.eq.3) then
445 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
446 else if (lang.eq.4) then
447 write (iout,'(a/)') " in overdamped mode"
449 write (iout,'(//a,i5)')
450 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
453 write (iout,'(a60,f10.5)') "Temperature:",t_bath
454 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
455 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
456 write (iout,'(a60,f10.5)')
457 & "Scaling factor of the friction forces:",scal_fric
458 if (surfarea) write (iout,'(2a,i10,a)')
459 & "Friction coefficients will be scaled by solvent-accessible",
460 & " surface area every",reset_fricmat," steps."
462 c Calculate friction coefficients and bounds of stochastic forces
463 eta=6*pi*cPoise*etawat
464 if(me.eq.king.or..not.out1file)
465 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
467 gamp=scal_fric*(pstok+rwat)*eta
468 stdfp=dsqrt(2*Rb*t_bath/d_time)
470 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
471 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
473 if(me.eq.king.or..not.out1file)then
474 write (iout,'(/2a/)')
475 & "Radii of site types and friction coefficients and std's of",
476 & " stochastic forces of fully exposed sites"
477 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
479 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
480 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
484 if(me.eq.king.or..not.out1file)then
485 write (iout,'(a)') "Berendsen bath calculation"
486 write (iout,'(a60,f10.5)') "Temperature:",t_bath
487 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
489 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
490 & count_reset_moment," steps"
492 & write (iout,'(a,i10,a)')
493 & "Velocities will be reset at random every",count_reset_vel,
496 else if (tnp .or. tnp1 .or. tnh) then
497 if (tnp .or. tnp1) then
498 write (iout,'(a)') "Nose-Poincare bath calculation"
499 if (tnp) write (iout,'(a)')
500 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
501 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
503 write (iout,'(a)') "Nose-Hoover bath calculation"
504 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
514 WDTI(i) = 1.0*d_time/nresn
521 write (iout,'(a)') "NVT-XI-RESPA algorithm"
523 write (iout,'(a)') "NVT-XO-RESPA algorithm"
526 WDTIi(i) = 1.0*d_time/nresn/ntime_split
534 write (iout,'(a60,f10.5)') "Temperature:",t_bath
535 write (iout,'(a60,f10.5)') "Q =",Q_np
537 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
538 & count_reset_moment," steps"
540 & write (iout,'(a,i10,a)')
541 & "Velocities will be reset at random every",count_reset_vel,
544 else if (hmc.gt.0) then
545 write (iout,'(a)') "Hybrid Monte Carlo calculation"
546 write (iout,'(a60,f10.5)') "Temperature:",t_bath
547 write (iout,'(a60,i10)')
548 & "Number of MD steps between Metropolis tests:",hmc
551 if(me.eq.king.or..not.out1file)
552 & write (iout,'(a31)') "Microcanonical mode calculation"
554 if(me.eq.king.or..not.out1file)then
555 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
557 write(iout,*) "MD running with constraints."
558 write(iout,*) "Equilibration time ", eq_time, " mtus."
559 write(iout,*) "Constraining ", nfrag," fragments."
560 write(iout,*) "Length of each fragment, weight and q0:"
562 write (iout,*) "Set of restraints #",iset
564 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
565 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
567 write(iout,*) "constraints between ", npair, "fragments."
568 write(iout,*) "constraint pairs, weights and q0:"
570 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
571 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
573 write(iout,*) "angle constraints within ", nfrag_back,
574 & "backbone fragments."
575 write(iout,*) "fragment, weights:"
577 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
578 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
579 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
582 iset=mod(kolor,nset)+1
585 if(me.eq.king.or..not.out1file)
586 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
589 c------------------------------------------------------------------------------
592 C Read molecular data.
594 implicit real*8 (a-h,o-z)
600 include 'COMMON.IOUNITS'
603 include 'COMMON.INTERACT'
604 include 'COMMON.LOCAL'
605 include 'COMMON.NAMES'
606 include 'COMMON.CHAIN'
607 include 'COMMON.FFIELD'
608 include 'COMMON.SBRIDGE'
609 include 'COMMON.HEADER'
610 include 'COMMON.CONTROL'
611 include 'COMMON.DBASE'
612 include 'COMMON.THREAD'
613 include 'COMMON.CONTACTS'
614 include 'COMMON.TORCNSTR'
615 include 'COMMON.TIME1'
616 include 'COMMON.BOUNDS'
618 include 'COMMON.REMD'
619 include 'COMMON.SETUP'
620 character*4 sequence(maxres)
622 double precision x(maxvar)
623 character*256 pdbfile
624 character*320 weightcard
625 character*80 weightcard_t,ucase
626 dimension itype_pdb(maxres)
627 common /pizda/ itype_pdb
628 logical seq_comp,fail
629 double precision energia(0:n_ene)
636 C Read weights of the subsequent energy terms.
649 if(me.eq.king.or..not.out1file) then
650 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
651 write (iout,*) 'Current weights for processor ',
652 & me,' set ',i2set(me)
656 call card_concat(weightcard)
657 call reada(weightcard,'WLONG',wlong,1.0D0)
658 call reada(weightcard,'WSC',wsc,wlong)
659 call reada(weightcard,'WSCP',wscp,wlong)
660 call reada(weightcard,'WELEC',welec,1.0D0)
661 call reada(weightcard,'WVDWPP',wvdwpp,welec)
662 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
663 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
664 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
665 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
666 call reada(weightcard,'WTURN3',wturn3,1.0D0)
667 call reada(weightcard,'WTURN4',wturn4,1.0D0)
668 call reada(weightcard,'WTURN6',wturn6,1.0D0)
669 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
670 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
671 call reada(weightcard,'WBOND',wbond,1.0D0)
672 call reada(weightcard,'WTOR',wtor,1.0D0)
673 call reada(weightcard,'WTORD',wtor_d,1.0D0)
674 call reada(weightcard,'WANG',wang,1.0D0)
675 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
676 call reada(weightcard,'SCAL14',scal14,0.4D0)
677 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
678 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
679 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
680 call reada(weightcard,'TEMP0',temp0,300.0d0)
681 if (index(weightcard,'SOFT').gt.0) ipot=6
682 C 12/1/95 Added weight for the multi-body term WCORR
683 call reada(weightcard,'WCORRH',wcorr,1.0D0)
684 if (wcorr4.gt.0.0d0) wcorr=wcorr4
692 hweights(i,7)=wel_loc
695 hweights(i,10)=wturn6
697 hweights(i,12)=wscloc
699 hweights(i,14)=wtor_d
700 hweights(i,15)=wstrain
701 hweights(i,16)=wvdwpp
703 hweights(i,18)=scal14
704 hweights(i,21)=wsccor
709 weights(i)=hweights(i2set(me),i)
733 call card_concat(weightcard)
734 call reada(weightcard,'WLONG',wlong,1.0D0)
735 call reada(weightcard,'WSC',wsc,wlong)
736 call reada(weightcard,'WSCP',wscp,wlong)
737 call reada(weightcard,'WELEC',welec,1.0D0)
738 call reada(weightcard,'WVDWPP',wvdwpp,welec)
739 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
740 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
741 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
742 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
743 call reada(weightcard,'WTURN3',wturn3,1.0D0)
744 call reada(weightcard,'WTURN4',wturn4,1.0D0)
745 call reada(weightcard,'WTURN6',wturn6,1.0D0)
746 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
747 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
748 call reada(weightcard,'WBOND',wbond,1.0D0)
749 call reada(weightcard,'WTOR',wtor,1.0D0)
750 call reada(weightcard,'WTORD',wtor_d,1.0D0)
751 call reada(weightcard,'WANG',wang,1.0D0)
752 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
753 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
754 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
755 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
756 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
757 call reada(weightcard,'SCAL14',scal14,0.4D0)
758 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
759 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
760 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
761 call reada(weightcard,'TEMP0',temp0,300.0d0)
762 if (index(weightcard,'SOFT').gt.0) ipot=6
763 C 12/1/95 Added weight for the multi-body term WCORR
764 call reada(weightcard,'WCORRH',wcorr,1.0D0)
765 if (wcorr4.gt.0.0d0) wcorr=wcorr4
786 weights(25)=wdfa_dist
789 weights(28)=wdfa_beta
791 if(me.eq.king.or..not.out1file)
792 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
793 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
795 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
796 10 format (/'Energy-term weights (unscaled):'//
797 & 'WSCC= ',f10.6,' (SC-SC)'/
798 & 'WSCP= ',f10.6,' (SC-p)'/
799 & 'WELEC= ',f10.6,' (p-p electr)'/
800 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
801 & 'WBOND= ',f10.6,' (stretching)'/
802 & 'WANG= ',f10.6,' (bending)'/
803 & 'WSCLOC= ',f10.6,' (SC local)'/
804 & 'WTOR= ',f10.6,' (torsional)'/
805 & 'WTORD= ',f10.6,' (double torsional)'/
806 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
807 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
808 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
809 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
810 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
811 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
812 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
813 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
814 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
815 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
816 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
817 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
818 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
819 if(me.eq.king.or..not.out1file)then
820 if (wcorr4.gt.0.0d0) then
821 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
822 & 'between contact pairs of peptide groups'
823 write (iout,'(2(a,f5.3/))')
824 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
825 & 'Range of quenching the correlation terms:',2*delt_corr
826 else if (wcorr.gt.0.0d0) then
827 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
828 & 'between contact pairs of peptide groups'
830 write (iout,'(a,f8.3)')
831 & 'Scaling factor of 1,4 SC-p interactions:',scal14
832 write (iout,'(a,f8.3)')
833 & 'General scaling factor of SC-p interactions:',scalscp
835 r0_corr=cutoff_corr-delt_corr
837 aad(i,1)=scalscp*aad(i,1)
838 aad(i,2)=scalscp*aad(i,2)
839 bad(i,1)=scalscp*bad(i,1)
840 bad(i,2)=scalscp*bad(i,2)
842 call rescale_weights(t_bath)
843 if(me.eq.king.or..not.out1file)
844 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
845 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
847 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
848 22 format (/'Energy-term weights (scaled):'//
849 & 'WSCC= ',f10.6,' (SC-SC)'/
850 & 'WSCP= ',f10.6,' (SC-p)'/
851 & 'WELEC= ',f10.6,' (p-p electr)'/
852 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
853 & 'WBOND= ',f10.6,' (stretching)'/
854 & 'WANG= ',f10.6,' (bending)'/
855 & 'WSCLOC= ',f10.6,' (SC local)'/
856 & 'WTOR= ',f10.6,' (torsional)'/
857 & 'WTORD= ',f10.6,' (double torsional)'/
858 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
859 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
860 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
861 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
862 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
863 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
864 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
865 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
866 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
867 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
868 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
869 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
870 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
871 if(me.eq.king.or..not.out1file)
872 & write (iout,*) "Reference temperature for weights calculation:",
874 call reada(weightcard,"D0CM",d0cm,3.78d0)
875 call reada(weightcard,"AKCM",akcm,15.1d0)
876 call reada(weightcard,"AKTH",akth,11.0d0)
877 call reada(weightcard,"AKCT",akct,12.0d0)
878 call reada(weightcard,"V1SS",v1ss,-1.08d0)
879 call reada(weightcard,"V2SS",v2ss,7.61d0)
880 call reada(weightcard,"V3SS",v3ss,13.7d0)
881 call reada(weightcard,"EBR",ebr,-5.50D0)
882 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
884 dyn_ss_mask(i)=.false.
888 dyn_ssbond_ij(i,j)=1.0d300
891 call reada(weightcard,"HT",Ht,0.0D0)
893 ss_depth=ebr/wsc-0.25*eps(1,1)
894 Ht=Ht/wsc-0.25*eps(1,1)
895 akcm=akcm*wstrain/wsc
896 akth=akth*wstrain/wsc
897 akct=akct*wstrain/wsc
898 v1ss=v1ss*wstrain/wsc
899 v2ss=v2ss*wstrain/wsc
900 v3ss=v3ss*wstrain/wsc
902 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
905 if(me.eq.king.or..not.out1file) then
906 write (iout,*) "Parameters of the SS-bond potential:"
907 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
909 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
910 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
911 write (iout,*)" HT",Ht
912 print *,'indpdb=',indpdb,' pdbref=',pdbref
914 if (indpdb.gt.0 .or. pdbref) then
915 read(inp,'(a)') pdbfile
916 if(me.eq.king.or..not.out1file)
917 & write (iout,'(2a)') 'PDB data will be read from file ',
918 & pdbfile(:ilen(pdbfile))
919 open(ipdbin,file=pdbfile,status='old',err=33)
921 33 write (iout,'(a)') 'Error opening PDB file.'
924 c print *,'Begin reading pdb data'
931 c print *,'Finished reading pdb data'
932 if(me.eq.king.or..not.out1file)
933 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
934 & ' nstart_sup=',nstart_sup
936 itype_pdb(i)=itype(i)
940 nct=nstart_sup+nsup-1
941 call contact(.false.,ncont_ref,icont_ref,co)
944 C Following 2 lines for diagnostics; comment out if not needed
945 write (iout,*) "Before sideadd"
947 if(me.eq.king.or..not.out1file)
948 & write(iout,*)'Adding sidechains'
955 do while (fail.and.nsi.le.maxsi)
956 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
959 if(fail) write(iout,*)'Adding sidechain failed for res ',
960 & i,' after ',nsi,' trials'
963 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
966 C Following 2 lines for diagnostics; comment out if not needed
967 c write (iout,*) "After sideadd"
970 if (indpdb.eq.0) then
971 C Read sequence if not taken from the pdb file.
973 c print *,'nres=',nres
974 if (iscode.gt.0) then
975 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
977 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
979 C Convert sequence to numeric code
981 itype(i)=rescode(i,sequence(i),iscode)
983 C Assign initial virtual bond lengths
989 vbld(i+nres)=dsc(itype(i))
990 vbld_inv(i+nres)=dsc_inv(itype(i))
991 c write (iout,*) "i",i," itype",itype(i),
992 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
996 c print '(20i4)',(itype(i),i=1,nres)
999 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1001 if (itype(i).eq.21) then
1005 else if (itype(i+1).ne.20) then
1007 else if (itype(i).ne.20) then
1014 if(me.eq.king.or..not.out1file)then
1015 write (iout,*) "ITEL"
1017 write (iout,*) i,itype(i),itel(i)
1019 print *,'Call Read_Bridge.'
1022 C 8/13/98 Set limits to generating the dihedral angles
1027 read (inp,*) ndih_constr
1028 if (ndih_constr.gt.0) then
1030 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1031 if(me.eq.king.or..not.out1file)then
1033 & 'There are',ndih_constr,' constraints on phi angles.'
1035 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1039 phi0(i)=deg2rad*phi0(i)
1040 drange(i)=deg2rad*drange(i)
1042 if(me.eq.king.or..not.out1file)
1043 & write (iout,*) 'FTORS',ftors
1046 phibound(1,ii) = phi0(i)-drange(i)
1047 phibound(2,ii) = phi0(i)+drange(i)
1052 if (me.eq.king) then
1054 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1056 write (iout,'(a3,i5,2f10.1)')
1057 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1063 cd print *,'NNT=',NNT,' NCT=',NCT
1064 if (itype(1).eq.21) nnt=2
1065 if (itype(nres).eq.21) nct=nct-1
1067 C Bartek:READ init_vars
1068 C Initialize variables!
1069 C Juyong:READ read_info
1070 C READ fragment information!!
1071 C both routines should be in dfa.F file!!
1073 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1074 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1076 print*, 'init_dfa_vars finished!'
1078 print*, 'read_dfa_info finished!'
1082 if(me.eq.king.or..not.out1file)
1083 & write (iout,'(a,i3)') 'nsup=',nsup
1085 if (nsup.le.(nct-nnt+1)) then
1086 do i=0,nct-nnt+1-nsup
1087 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1093 & 'Error - sequences to be superposed do not match.'
1096 do i=0,nsup-(nct-nnt+1)
1097 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1099 nstart_sup=nstart_sup+i
1105 & 'Error - sequences to be superposed do not match.'
1108 if (nsup.eq.0) nsup=nct-nnt
1109 if (nstart_sup.eq.0) nstart_sup=nnt
1110 if (nstart_seq.eq.0) nstart_seq=nnt
1111 if(me.eq.king.or..not.out1file)
1112 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1113 & ' nstart_seq=',nstart_seq
1115 c--- Zscore rms -------
1116 if (nz_start.eq.0) nz_start=nnt
1117 if (nz_end.eq.0 .and. nsup.gt.0) then
1119 else if (nz_end.eq.0) then
1122 if(me.eq.king.or..not.out1file)then
1123 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1124 write (iout,*) 'IZ_SC=',iz_sc
1126 c----------------------
1129 if (.not.pdbref) then
1130 call read_angles(inp,*38)
1132 38 write (iout,'(a)') 'Error reading reference structure.'
1134 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1135 stop 'Error reading reference structure'
1139 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1148 call contact(.true.,ncont_ref,icont_ref,co)
1150 if(me.eq.king.or..not.out1file)
1151 & write (iout,*) 'Contact order:',co
1153 if(me.eq.king.or..not.out1file)
1154 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1157 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1159 if(me.eq.king.or..not.out1file)
1160 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1161 & icont_ref(1,i),' ',
1162 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1166 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1167 if (constr_dist.gt.0) then
1168 call read_dist_constr
1172 if (constr_homology.gt.0) then
1173 call read_constr_homology
1174 if (indpdb.gt.0 .or. pdbref) then
1177 c(j,i)=crefjlee(j,i)
1178 cref(j,i)=crefjlee(j,i)
1187 if (nhpb.gt.0) call hpb_partition
1188 c write (iout,*) "After read_dist_constr nhpb",nhpb
1190 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1191 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1192 & modecalc.ne.10) then
1193 C If input structure hasn't been supplied from the PDB file read or generate
1195 if (iranconf.eq.0 .and. .not. extconf) then
1196 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1197 & write (iout,'(a)') 'Initial geometry will be read in.'
1199 read(inp,'(8f10.5)',end=36,err=36)
1200 & ((c(l,k),l=1,3),k=1,nres),
1201 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1202 call int_from_cart1(.false.)
1205 dc(j,i)=c(j,i+1)-c(j,i)
1206 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1210 if (itype(i).ne.10) then
1212 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1213 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1219 call read_angles(inp,*36)
1222 36 write (iout,'(a)') 'Error reading angle file.'
1224 call mpi_finalize( MPI_COMM_WORLD,IERR )
1226 stop 'Error reading angle file.'
1228 else if (extconf) then
1229 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1230 & write (iout,'(a)') 'Extended chain initial geometry.'
1232 theta(i)=90d0*deg2rad
1235 phi(i)=180d0*deg2rad
1238 alph(i)=110d0*deg2rad
1241 omeg(i)=-120d0*deg2rad
1244 if(me.eq.king.or..not.out1file)
1245 & write (iout,'(a)') 'Random-generated initial geometry.'
1249 if (me.eq.king .or. fg_rank.eq.0 .and. (
1250 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1254 call gen_rand_conf(itmp,*30)
1256 30 write (iout,*) 'Failed to generate random conformation',
1257 & ', itrial=',itrial
1258 write (*,*) 'Processor:',me,
1259 & ' Failed to generate random conformation',
1269 write (iout,'(a,i3,a)') 'Processor:',me,
1270 & ' error in generating random conformation.'
1271 write (*,'(a,i3,a)') 'Processor:',me,
1272 & ' error in generating random conformation.'
1275 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1282 elseif (modecalc.eq.4) then
1283 read (inp,'(a)') intinname
1284 open (intin,file=intinname,status='old',err=333)
1285 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1286 & write (iout,'(a)') 'intinname',intinname
1287 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1289 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1291 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1293 stop 'Error opening angle file.'
1297 C Generate distance constraints, if the PDB structure is to be regularized.
1298 if (nthread.gt.0) then
1299 call read_threadbase
1302 if (me.eq.king .or. .not. out1file)
1304 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1305 write (iout,'(/a,i3,a)')
1306 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1307 write (iout,'(20i4)') (iss(i),i=1,ns)
1309 write(iout,*)"Running with dynamic disulfide-bond formation"
1311 write (iout,'(/a/)') 'Pre-formed links are:'
1317 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1318 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1324 if (ns.gt.0.and.dyn_ss) then
1328 forcon(i-nss)=forcon(i)
1335 dyn_ss_mask(iss(i))=.true.
1338 if (i2ndstr.gt.0) call secstrp2dihc
1339 c call geom_to_var(nvar,x)
1340 c call etotal(energia(0))
1341 c call enerprint(energia(0))
1342 c call briefout(0,etot)
1344 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1345 cd write (iout,'(a)') 'Variable list:'
1346 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1348 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1349 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1350 & 'Processor',myrank,': end reading molecular data.'
1354 c--------------------------------------------------------------------------
1355 logical function seq_comp(itypea,itypeb,length)
1357 integer length,itypea(length),itypeb(length)
1360 if (itypea(i).ne.itypeb(i)) then
1368 c-----------------------------------------------------------------------------
1369 subroutine read_bridge
1370 C Read information about disulfide bridges.
1371 implicit real*8 (a-h,o-z)
1372 include 'DIMENSIONS'
1376 include 'COMMON.IOUNITS'
1377 include 'COMMON.GEO'
1378 include 'COMMON.VAR'
1379 include 'COMMON.INTERACT'
1380 include 'COMMON.LOCAL'
1381 include 'COMMON.NAMES'
1382 include 'COMMON.CHAIN'
1383 include 'COMMON.FFIELD'
1384 include 'COMMON.SBRIDGE'
1385 include 'COMMON.HEADER'
1386 include 'COMMON.CONTROL'
1387 include 'COMMON.DBASE'
1388 include 'COMMON.THREAD'
1389 include 'COMMON.TIME1'
1390 include 'COMMON.SETUP'
1391 C Read bridging residues.
1392 read (inp,*) ns,(iss(i),i=1,ns)
1394 if(me.eq.king.or..not.out1file)
1395 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1396 C Check whether the specified bridging residues are cystines.
1398 if (itype(iss(i)).ne.1) then
1399 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1400 & 'Do you REALLY think that the residue ',
1401 & restyp(itype(iss(i))),i,
1402 & ' can form a disulfide bridge?!!!'
1403 write (*,'(2a,i3,a)')
1404 & 'Do you REALLY think that the residue ',
1405 & restyp(itype(iss(i))),i,
1406 & ' can form a disulfide bridge?!!!'
1408 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1413 C Read preformed bridges.
1415 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1417 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1420 C Check if the residues involved in bridges are in the specified list of
1421 C bridging residues.
1424 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1425 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1426 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1427 & ' contains residues present in other pairs.'
1428 write (*,'(a,i3,a)') 'Disulfide pair',i,
1429 & ' contains residues present in other pairs.'
1431 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1437 if (ihpb(i).eq.iss(j)) goto 10
1439 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1442 if (jhpb(i).eq.iss(j)) goto 20
1444 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1450 ihpb(i)=ihpb(i)+nres
1451 jhpb(i)=jhpb(i)+nres
1457 c----------------------------------------------------------------------------
1458 subroutine read_x(kanal,*)
1459 implicit real*8 (a-h,o-z)
1460 include 'DIMENSIONS'
1461 include 'COMMON.GEO'
1462 include 'COMMON.VAR'
1463 include 'COMMON.CHAIN'
1464 include 'COMMON.IOUNITS'
1465 include 'COMMON.CONTROL'
1466 include 'COMMON.LOCAL'
1467 include 'COMMON.INTERACT'
1468 c Read coordinates from input
1470 read(kanal,'(8f10.5)',end=10,err=10)
1471 & ((c(l,k),l=1,3),k=1,nres),
1472 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1475 c(j,2*nres)=c(j,nres)
1477 call int_from_cart1(.false.)
1480 dc(j,i)=c(j,i+1)-c(j,i)
1481 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1485 if (itype(i).ne.10) then
1487 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1488 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1496 c----------------------------------------------------------------------------
1497 subroutine read_threadbase
1498 implicit real*8 (a-h,o-z)
1499 include 'DIMENSIONS'
1500 include 'COMMON.IOUNITS'
1501 include 'COMMON.GEO'
1502 include 'COMMON.VAR'
1503 include 'COMMON.INTERACT'
1504 include 'COMMON.LOCAL'
1505 include 'COMMON.NAMES'
1506 include 'COMMON.CHAIN'
1507 include 'COMMON.FFIELD'
1508 include 'COMMON.SBRIDGE'
1509 include 'COMMON.HEADER'
1510 include 'COMMON.CONTROL'
1511 include 'COMMON.DBASE'
1512 include 'COMMON.THREAD'
1513 include 'COMMON.TIME1'
1514 C Read pattern database for threading.
1515 read (icbase,*) nseq
1517 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1518 & nres_base(2,i),nres_base(3,i)
1519 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1521 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1522 c & nres_base(2,i),nres_base(3,i)
1523 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1527 if (weidis.eq.0.0D0) weidis=0.1D0
1536 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1537 write (iout,'(a,i5)') 'nexcl: ',nexcl
1538 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1541 c------------------------------------------------------------------------------
1542 subroutine setup_var
1543 implicit real*8 (a-h,o-z)
1544 include 'DIMENSIONS'
1545 include 'COMMON.IOUNITS'
1546 include 'COMMON.GEO'
1547 include 'COMMON.VAR'
1548 include 'COMMON.INTERACT'
1549 include 'COMMON.LOCAL'
1550 include 'COMMON.NAMES'
1551 include 'COMMON.CHAIN'
1552 include 'COMMON.FFIELD'
1553 include 'COMMON.SBRIDGE'
1554 include 'COMMON.HEADER'
1555 include 'COMMON.CONTROL'
1556 include 'COMMON.DBASE'
1557 include 'COMMON.THREAD'
1558 include 'COMMON.TIME1'
1559 C Set up variable list.
1565 if (itype(i).ne.10) then
1567 ialph(i,1)=nvar+nside
1571 if (indphi.gt.0) then
1573 else if (indback.gt.0) then
1578 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1581 c----------------------------------------------------------------------------
1582 subroutine gen_dist_constr
1583 C Generate CA distance constraints.
1584 implicit real*8 (a-h,o-z)
1585 include 'DIMENSIONS'
1586 include 'COMMON.IOUNITS'
1587 include 'COMMON.GEO'
1588 include 'COMMON.VAR'
1589 include 'COMMON.INTERACT'
1590 include 'COMMON.LOCAL'
1591 include 'COMMON.NAMES'
1592 include 'COMMON.CHAIN'
1593 include 'COMMON.FFIELD'
1594 include 'COMMON.SBRIDGE'
1595 include 'COMMON.HEADER'
1596 include 'COMMON.CONTROL'
1597 include 'COMMON.DBASE'
1598 include 'COMMON.THREAD'
1599 include 'COMMON.TIME1'
1600 dimension itype_pdb(maxres)
1601 common /pizda/ itype_pdb
1603 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1604 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1605 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1607 do i=nstart_sup,nstart_sup+nsup-1
1608 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1609 cd & ' seq_pdb', restyp(itype_pdb(i))
1610 do j=i+2,nstart_sup+nsup-1
1612 ihpb(nhpb)=i+nstart_seq-nstart_sup
1613 jhpb(nhpb)=j+nstart_seq-nstart_sup
1615 dhpb(nhpb)=dist(i,j)
1618 cd write (iout,'(a)') 'Distance constraints:'
1623 cd if (ii.gt.nres) then
1628 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1629 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1630 cd & dhpb(i),forcon(i)
1634 c----------------------------------------------------------------------------
1636 implicit real*8 (a-h,o-z)
1637 include 'DIMENSIONS'
1638 include 'COMMON.MAP'
1639 include 'COMMON.IOUNITS'
1640 character*3 angid(4) /'THE','PHI','ALP','OME'/
1641 character*80 mapcard,ucase
1643 read (inp,'(a)') mapcard
1644 mapcard=ucase(mapcard)
1645 if (index(mapcard,'PHI').gt.0) then
1647 else if (index(mapcard,'THE').gt.0) then
1649 else if (index(mapcard,'ALP').gt.0) then
1651 else if (index(mapcard,'OME').gt.0) then
1654 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1655 stop 'Error - illegal variable spec in MAP card.'
1657 call readi (mapcard,'RES1',res1(imap),0)
1658 call readi (mapcard,'RES2',res2(imap),0)
1659 if (res1(imap).eq.0) then
1660 res1(imap)=res2(imap)
1661 else if (res2(imap).eq.0) then
1662 res2(imap)=res1(imap)
1664 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1666 & 'Error - illegal definition of variable group in MAP.'
1667 stop 'Error - illegal definition of variable group in MAP.'
1669 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1670 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1671 call readi(mapcard,'NSTEP',nstep(imap),0)
1672 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1674 & 'Illegal boundary and/or step size specification in MAP.'
1675 stop 'Illegal boundary and/or step size specification in MAP.'
1680 c----------------------------------------------------------------------------
1681 csa subroutine csaread
1682 csa implicit real*8 (a-h,o-z)
1683 csa include 'DIMENSIONS'
1684 csa include 'COMMON.IOUNITS'
1685 csa include 'COMMON.GEO'
1686 csa include 'COMMON.CSA'
1687 csa include 'COMMON.BANK'
1688 csa include 'COMMON.CONTROL'
1689 csa character*80 ucase
1690 csa character*620 mcmcard
1691 csa call card_concat(mcmcard)
1693 csa call readi(mcmcard,'NCONF',nconf,50)
1694 csa call readi(mcmcard,'NADD',nadd,0)
1695 csa call readi(mcmcard,'JSTART',jstart,1)
1696 csa call readi(mcmcard,'JEND',jend,1)
1697 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1698 csa call readi(mcmcard,'N0',n0,1)
1699 csa call readi(mcmcard,'N1',n1,6)
1700 csa call readi(mcmcard,'N2',n2,4)
1701 csa call readi(mcmcard,'N3',n3,0)
1702 csa call readi(mcmcard,'N4',n4,0)
1703 csa call readi(mcmcard,'N5',n5,0)
1704 csa call readi(mcmcard,'N6',n6,10)
1705 csa call readi(mcmcard,'N7',n7,0)
1706 csa call readi(mcmcard,'N8',n8,0)
1707 csa call readi(mcmcard,'N9',n9,0)
1708 csa call readi(mcmcard,'N14',n14,0)
1709 csa call readi(mcmcard,'N15',n15,0)
1710 csa call readi(mcmcard,'N16',n16,0)
1711 csa call readi(mcmcard,'N17',n17,0)
1712 csa call readi(mcmcard,'N18',n18,0)
1714 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1716 csa call readi(mcmcard,'NDIFF',ndiff,2)
1717 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1718 csa call readi(mcmcard,'IS1',is1,1)
1719 csa call readi(mcmcard,'IS2',is2,8)
1720 csa call readi(mcmcard,'NRAN0',nran0,4)
1721 csa call readi(mcmcard,'NRAN1',nran1,2)
1722 csa call readi(mcmcard,'IRR',irr,1)
1723 csa call readi(mcmcard,'NSEED',nseed,20)
1724 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1725 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1726 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1727 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1728 csa call readi(mcmcard,'ICMAX',icmax,3)
1729 csa call readi(mcmcard,'IRESTART',irestart,0)
1730 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1733 csa call reada(mcmcard,'DELE',dele,20.0d0)
1734 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1735 csa call readi(mcmcard,'IREF',iref,0)
1736 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1737 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1738 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1739 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1740 csa write (iout,*) "NCONF_IN",nconf_in
1743 c----------------------------------------------------------------------------
1744 cfmc subroutine mcmfread
1745 cfmc implicit real*8 (a-h,o-z)
1746 cfmc include 'DIMENSIONS'
1747 cfmc include 'COMMON.MCMF'
1748 cfmc include 'COMMON.IOUNITS'
1749 cfmc include 'COMMON.GEO'
1750 cfmc character*80 ucase
1751 cfmc character*620 mcmcard
1752 cfmc call card_concat(mcmcard)
1754 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1755 cfmc write(iout,*)'MAXRANT=',maxrant
1756 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1757 cfmc write(iout,*)'MAXFAM=',maxfam
1758 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1759 cfmc write(iout,*)'NNET1=',nnet1
1760 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1761 cfmc write(iout,*)'NNET2=',nnet2
1762 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1763 cfmc write(iout,*)'NNET3=',nnet3
1764 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1765 cfmc write(iout,*)'ILASTT=',ilastt
1766 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1767 cfmc write(iout,*)'MAXSTR=',maxstr
1768 cfmc maxstr_f=maxstr/maxfam
1769 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1770 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1771 cfmc write(iout,*)'NMCMF=',nmcmf
1772 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1773 cfmc write(iout,*)'IFOCUS=',ifocus
1774 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1775 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1776 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1777 cfmc write(iout,*)'INTPRT=',intprt
1778 cfmc call readi(mcmcard,'IPRT',iprt,100)
1779 cfmc write(iout,*)'IPRT=',iprt
1780 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1781 cfmc write(iout,*)'IMAXTR=',imaxtr
1782 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1783 cfmc write(iout,*)'MAXEVEN=',maxeven
1784 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1785 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1786 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1787 cfmc write(iout,*)'INIMIN=',inimin
1788 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1789 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1790 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1791 cfmc write(iout,*)'NTHREAD=',nthread
1792 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1793 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1794 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1795 cfmc write(iout,*)'MAXPERT=',maxpert
1796 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1797 cfmc write(iout,*)'IRMSD=',irmsd
1798 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1799 cfmc write(iout,*)'DENEMIN=',denemin
1800 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1801 cfmc write(iout,*)'RCUT1S=',rcut1s
1802 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1803 cfmc write(iout,*)'RCUT1E=',rcut1e
1804 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1805 cfmc write(iout,*)'RCUT2S=',rcut2s
1806 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1807 cfmc write(iout,*)'RCUT2E=',rcut2e
1808 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1809 cfmc write(iout,*)'DPERT1=',d_pert1
1810 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1811 cfmc write(iout,*)'DPERT1A=',d_pert1a
1812 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1813 cfmc write(iout,*)'DPERT2=',d_pert2
1814 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1815 cfmc write(iout,*)'DPERT2A=',d_pert2a
1816 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1817 cfmc write(iout,*)'DPERT2B=',d_pert2b
1818 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1819 cfmc write(iout,*)'DPERT2C=',d_pert2c
1820 cfmc d_pert1=deg2rad*d_pert1
1821 cfmc d_pert1a=deg2rad*d_pert1a
1822 cfmc d_pert2=deg2rad*d_pert2
1823 cfmc d_pert2a=deg2rad*d_pert2a
1824 cfmc d_pert2b=deg2rad*d_pert2b
1825 cfmc d_pert2c=deg2rad*d_pert2c
1826 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1827 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1828 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1829 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1830 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1831 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1832 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1833 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1834 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1835 cfmc write(iout,*)'RCUTINI=',rcutini
1836 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1837 cfmc write(iout,*)'GRAT=',grat
1838 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1839 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1843 c----------------------------------------------------------------------------
1845 implicit real*8 (a-h,o-z)
1846 include 'DIMENSIONS'
1847 include 'COMMON.MCM'
1848 include 'COMMON.MCE'
1849 include 'COMMON.IOUNITS'
1851 character*320 mcmcard
1852 call card_concat(mcmcard)
1853 call readi(mcmcard,'MAXACC',maxacc,100)
1854 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1855 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1856 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1857 call readi(mcmcard,'MAXREPM',maxrepm,200)
1858 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1859 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1860 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1861 call reada(mcmcard,'E_UP',e_up,5.0D0)
1862 call reada(mcmcard,'DELTE',delte,0.1D0)
1863 call readi(mcmcard,'NSWEEP',nsweep,5)
1864 call readi(mcmcard,'NSTEPH',nsteph,0)
1865 call readi(mcmcard,'NSTEPC',nstepc,0)
1866 call reada(mcmcard,'TMIN',tmin,298.0D0)
1867 call reada(mcmcard,'TMAX',tmax,298.0D0)
1868 call readi(mcmcard,'NWINDOW',nwindow,0)
1869 call readi(mcmcard,'PRINT_MC',print_mc,0)
1870 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1871 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1872 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1873 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1874 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1875 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1876 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1877 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1878 if (nwindow.gt.0) then
1879 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1881 winlen(i)=winend(i)-winstart(i)+1
1884 if (tmax.lt.tmin) tmax=tmin
1885 if (tmax.eq.tmin) then
1889 if (nstepc.gt.0 .and. nsteph.gt.0) then
1890 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1891 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1893 C Probabilities of different move types
1894 sumpro_type(0)=0.0D0
1895 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1896 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1897 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1898 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1899 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1900 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1901 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1903 print *,'i',i,' sumprotype',sumpro_type(i)
1904 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1905 print *,'i',i,' sumprotype',sumpro_type(i)
1909 c----------------------------------------------------------------------------
1910 subroutine read_minim
1911 implicit real*8 (a-h,o-z)
1912 include 'DIMENSIONS'
1913 include 'COMMON.MINIM'
1914 include 'COMMON.IOUNITS'
1916 character*320 minimcard
1917 call card_concat(minimcard)
1918 call readi(minimcard,'MAXMIN',maxmin,2000)
1919 call readi(minimcard,'MAXFUN',maxfun,5000)
1920 call readi(minimcard,'MINMIN',minmin,maxmin)
1921 call readi(minimcard,'MINFUN',minfun,maxmin)
1922 call reada(minimcard,'TOLF',tolf,1.0D-2)
1923 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1924 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1925 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1926 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1927 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1928 & 'Options in energy minimization:'
1929 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1930 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1931 & 'MinMin:',MinMin,' MinFun:',MinFun,
1932 & ' TolF:',TolF,' RTolF:',RTolF
1935 c----------------------------------------------------------------------------
1936 subroutine read_angles(kanal,*)
1937 implicit real*8 (a-h,o-z)
1938 include 'DIMENSIONS'
1939 include 'COMMON.GEO'
1940 include 'COMMON.VAR'
1941 include 'COMMON.CHAIN'
1942 include 'COMMON.IOUNITS'
1943 include 'COMMON.CONTROL'
1944 c Read angles from input
1946 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1947 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1948 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1949 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1952 c 9/7/01 avoid 180 deg valence angle
1953 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1955 theta(i)=deg2rad*theta(i)
1956 phi(i)=deg2rad*phi(i)
1957 alph(i)=deg2rad*alph(i)
1958 omeg(i)=deg2rad*omeg(i)
1963 c----------------------------------------------------------------------------
1964 subroutine reada(rekord,lancuch,wartosc,default)
1966 character*(*) rekord,lancuch
1967 double precision wartosc,default
1970 iread=index(rekord,lancuch)
1971 if (iread.eq.0) then
1975 iread=iread+ilen(lancuch)+1
1976 read (rekord(iread:),*,err=10,end=10) wartosc
1981 c----------------------------------------------------------------------------
1982 subroutine readi(rekord,lancuch,wartosc,default)
1984 character*(*) rekord,lancuch
1985 integer wartosc,default
1988 iread=index(rekord,lancuch)
1989 if (iread.eq.0) then
1993 iread=iread+ilen(lancuch)+1
1994 read (rekord(iread:),*,err=10,end=10) wartosc
1999 c----------------------------------------------------------------------------
2000 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2003 integer tablica(dim),default
2004 character*(*) rekord,lancuch
2011 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2012 if (iread.eq.0) return
2013 iread=iread+ilen(lancuch)+1
2014 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2017 c----------------------------------------------------------------------------
2018 subroutine multreada(rekord,lancuch,tablica,dim,default)
2021 double precision tablica(dim),default
2022 character*(*) rekord,lancuch
2029 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2030 if (iread.eq.0) return
2031 iread=iread+ilen(lancuch)+1
2032 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2035 c----------------------------------------------------------------------------
2036 subroutine openunits
2037 implicit real*8 (a-h,o-z)
2038 include 'DIMENSIONS'
2041 character*16 form,nodename
2044 include 'COMMON.SETUP'
2045 include 'COMMON.IOUNITS'
2047 include 'COMMON.CONTROL'
2048 integer lenpre,lenpot,ilen,lentmp
2050 character*3 out1file_text,ucase
2053 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2054 call getenv_loc("PREFIX",prefix)
2056 call getenv_loc("POT",pot)
2057 call getenv_loc("DIRTMP",tmpdir)
2058 call getenv_loc("CURDIR",curdir)
2059 call getenv_loc("OUT1FILE",out1file_text)
2060 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2061 out1file_text=ucase(out1file_text)
2062 if (out1file_text(1:1).eq."Y") then
2065 out1file=fg_rank.gt.0
2070 if (lentmp.gt.0) then
2071 write (*,'(80(1h!))')
2072 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2073 write (*,'(80(1h!))')
2074 write (*,*)"All output files will be on node /tmp directory."
2076 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2077 if (me.eq.king) then
2078 write (*,*) "The master node is ",nodename
2079 else if (fg_rank.eq.0) then
2080 write (*,*) "I am the CG slave node ",nodename
2082 write (*,*) "I am the FG slave node ",nodename
2085 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2086 lenpre = lentmp+lenpre+1
2088 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2089 C Get the names and open the input files
2090 #if defined(WINIFL) || defined(WINPGI)
2091 open(1,file=pref_orig(:ilen(pref_orig))//
2092 & '.inp',status='old',readonly,shared)
2093 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2094 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2095 C Get parameter filenames and open the parameter files.
2096 call getenv_loc('BONDPAR',bondname)
2097 open (ibond,file=bondname,status='old',readonly,shared)
2098 call getenv_loc('THETPAR',thetname)
2099 open (ithep,file=thetname,status='old',readonly,shared)
2101 call getenv_loc('THETPARPDB',thetname_pdb)
2102 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2104 call getenv_loc('ROTPAR',rotname)
2105 open (irotam,file=rotname,status='old',readonly,shared)
2107 call getenv_loc('ROTPARPDB',rotname_pdb)
2108 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2110 call getenv_loc('TORPAR',torname)
2111 open (itorp,file=torname,status='old',readonly,shared)
2112 call getenv_loc('TORDPAR',tordname)
2113 open (itordp,file=tordname,status='old',readonly,shared)
2114 call getenv_loc('FOURIER',fouriername)
2115 open (ifourier,file=fouriername,status='old',readonly,shared)
2116 call getenv_loc('ELEPAR',elename)
2117 open (ielep,file=elename,status='old',readonly,shared)
2118 call getenv_loc('SIDEPAR',sidename)
2119 open (isidep,file=sidename,status='old',readonly,shared)
2120 #elif (defined CRAY) || (defined AIX)
2121 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2123 c print *,"Processor",myrank," opened file 1"
2124 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2125 c print *,"Processor",myrank," opened file 9"
2126 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2127 C Get parameter filenames and open the parameter files.
2128 call getenv_loc('BONDPAR',bondname)
2129 open (ibond,file=bondname,status='old',action='read')
2130 c print *,"Processor",myrank," opened file IBOND"
2131 call getenv_loc('THETPAR',thetname)
2132 open (ithep,file=thetname,status='old',action='read')
2133 c print *,"Processor",myrank," opened file ITHEP"
2135 call getenv_loc('THETPARPDB',thetname_pdb)
2136 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2138 call getenv_loc('ROTPAR',rotname)
2139 open (irotam,file=rotname,status='old',action='read')
2140 c print *,"Processor",myrank," opened file IROTAM"
2142 call getenv_loc('ROTPARPDB',rotname_pdb)
2143 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2145 call getenv_loc('TORPAR',torname)
2146 open (itorp,file=torname,status='old',action='read')
2147 c print *,"Processor",myrank," opened file ITORP"
2148 call getenv_loc('TORDPAR',tordname)
2149 open (itordp,file=tordname,status='old',action='read')
2150 c print *,"Processor",myrank," opened file ITORDP"
2151 call getenv_loc('SCCORPAR',sccorname)
2152 open (isccor,file=sccorname,status='old',action='read')
2153 c print *,"Processor",myrank," opened file ISCCOR"
2154 call getenv_loc('FOURIER',fouriername)
2155 open (ifourier,file=fouriername,status='old',action='read')
2156 c print *,"Processor",myrank," opened file IFOURIER"
2157 call getenv_loc('ELEPAR',elename)
2158 open (ielep,file=elename,status='old',action='read')
2159 c print *,"Processor",myrank," opened file IELEP"
2160 call getenv_loc('SIDEPAR',sidename)
2161 open (isidep,file=sidename,status='old',action='read')
2162 c print *,"Processor",myrank," opened file ISIDEP"
2163 c print *,"Processor",myrank," opened parameter files"
2165 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2166 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2167 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2168 C Get parameter filenames and open the parameter files.
2169 call getenv_loc('BONDPAR',bondname)
2170 open (ibond,file=bondname,status='old')
2171 call getenv_loc('THETPAR',thetname)
2172 open (ithep,file=thetname,status='old')
2174 call getenv_loc('THETPARPDB',thetname_pdb)
2175 open (ithep_pdb,file=thetname_pdb,status='old')
2177 call getenv_loc('ROTPAR',rotname)
2178 open (irotam,file=rotname,status='old')
2180 call getenv_loc('ROTPARPDB',rotname_pdb)
2181 open (irotam_pdb,file=rotname_pdb,status='old')
2183 call getenv_loc('TORPAR',torname)
2184 open (itorp,file=torname,status='old')
2185 call getenv_loc('TORDPAR',tordname)
2186 open (itordp,file=tordname,status='old')
2187 call getenv_loc('SCCORPAR',sccorname)
2188 open (isccor,file=sccorname,status='old')
2189 call getenv_loc('FOURIER',fouriername)
2190 open (ifourier,file=fouriername,status='old')
2191 call getenv_loc('ELEPAR',elename)
2192 open (ielep,file=elename,status='old')
2193 call getenv_loc('SIDEPAR',sidename)
2194 open (isidep,file=sidename,status='old')
2196 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2198 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2199 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2200 C Get parameter filenames and open the parameter files.
2201 call getenv_loc('BONDPAR',bondname)
2202 open (ibond,file=bondname,status='old',action='read')
2203 call getenv_loc('THETPAR',thetname)
2204 open (ithep,file=thetname,status='old',action='read')
2206 call getenv_loc('THETPARPDB',thetname_pdb)
2207 print *,"thetname_pdb ",thetname_pdb
2208 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2209 print *,ithep_pdb," opened"
2211 call getenv_loc('ROTPAR',rotname)
2212 open (irotam,file=rotname,status='old',action='read')
2214 call getenv_loc('ROTPARPDB',rotname_pdb)
2215 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2217 call getenv_loc('TORPAR',torname)
2218 open (itorp,file=torname,status='old',action='read')
2219 call getenv_loc('TORDPAR',tordname)
2220 open (itordp,file=tordname,status='old',action='read')
2221 call getenv_loc('SCCORPAR',sccorname)
2222 open (isccor,file=sccorname,status='old',action='read')
2223 call getenv_loc('FOURIER',fouriername)
2224 open (ifourier,file=fouriername,status='old',action='read')
2225 call getenv_loc('ELEPAR',elename)
2226 open (ielep,file=elename,status='old',action='read')
2227 call getenv_loc('SIDEPAR',sidename)
2228 open (isidep,file=sidename,status='old',action='read')
2232 C 8/9/01 In the newest version SCp interaction constants are read from a file
2233 C Use -DOLDSCP to use hard-coded constants instead.
2235 call getenv_loc('SCPPAR',scpname)
2236 #if defined(WINIFL) || defined(WINPGI)
2237 open (iscpp,file=scpname,status='old',readonly,shared)
2238 #elif (defined CRAY) || (defined AIX)
2239 open (iscpp,file=scpname,status='old',action='read')
2241 open (iscpp,file=scpname,status='old')
2243 open (iscpp,file=scpname,status='old',action='read')
2246 call getenv_loc('PATTERN',patname)
2247 #if defined(WINIFL) || defined(WINPGI)
2248 open (icbase,file=patname,status='old',readonly,shared)
2249 #elif (defined CRAY) || (defined AIX)
2250 open (icbase,file=patname,status='old',action='read')
2252 open (icbase,file=patname,status='old')
2254 open (icbase,file=patname,status='old',action='read')
2257 C Open output file only for CG processes
2258 c print *,"Processor",myrank," fg_rank",fg_rank
2259 if (fg_rank.eq.0) then
2261 if (nodes.eq.1) then
2264 npos = dlog10(dfloat(nodes-1))+1
2266 if (npos.lt.3) npos=3
2267 write (liczba,'(i1)') npos
2268 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2270 write (liczba,form) me
2271 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2272 & liczba(:ilen(liczba))
2273 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2275 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2277 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2278 & liczba(:ilen(liczba))//'.mol2'
2279 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2280 & liczba(:ilen(liczba))//'.stat'
2282 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2283 & //liczba(:ilen(liczba))//'.stat')
2284 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2287 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2288 & liczba(:ilen(liczba))//'.const'
2293 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2294 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2295 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2296 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2297 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2299 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2301 rest2name=prefix(:ilen(prefix))//'.rst'
2303 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2306 #if defined(AIX) || defined(PGI)
2307 if (me.eq.king .or. .not. out1file)
2308 & open(iout,file=outname,status='unknown')
2311 if (fg_rank.gt.0) then
2312 write (liczba,'(i3.3)') myrank/nfgtasks
2313 write (ll,'(bz,i3.3)') fg_rank
2314 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2320 open(igeom,file=intname,status='unknown',position='append')
2321 open(ipdb,file=pdbname,status='unknown')
2322 open(imol2,file=mol2name,status='unknown')
2323 open(istat,file=statname,status='unknown',position='append')
2325 c1out open(iout,file=outname,status='unknown')
2328 if (me.eq.king .or. .not.out1file)
2329 & open(iout,file=outname,status='unknown')
2332 if (fg_rank.gt.0) then
2333 print "Processor",fg_rank," opening output file"
2334 write (liczba,'(i3.3)') myrank/nfgtasks
2335 write (ll,'(bz,i3.3)') fg_rank
2336 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2342 open(igeom,file=intname,status='unknown',access='append')
2343 open(ipdb,file=pdbname,status='unknown')
2344 open(imol2,file=mol2name,status='unknown')
2345 open(istat,file=statname,status='unknown',access='append')
2347 c1out open(iout,file=outname,status='unknown')
2350 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2351 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2352 csa csa_history=prefix(:lenpre)//'.CSA.history'
2353 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2354 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2355 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2356 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2357 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2358 csa csa_int=prefix(:lenpre)//'.int'
2359 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2360 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2361 csa csa_in=prefix(:lenpre)//'.CSA.in'
2362 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2365 write (iout,'(80(1h-))')
2366 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2367 write (iout,'(80(1h-))')
2368 write (iout,*) "Input file : ",
2369 & pref_orig(:ilen(pref_orig))//'.inp'
2370 write (iout,*) "Output file : ",
2371 & outname(:ilen(outname))
2373 write (iout,*) "Sidechain potential file : ",
2374 & sidename(:ilen(sidename))
2376 write (iout,*) "SCp potential file : ",
2377 & scpname(:ilen(scpname))
2379 write (iout,*) "Electrostatic potential file : ",
2380 & elename(:ilen(elename))
2381 write (iout,*) "Cumulant coefficient file : ",
2382 & fouriername(:ilen(fouriername))
2383 write (iout,*) "Torsional parameter file : ",
2384 & torname(:ilen(torname))
2385 write (iout,*) "Double torsional parameter file : ",
2386 & tordname(:ilen(tordname))
2387 write (iout,*) "SCCOR parameter file : ",
2388 & sccorname(:ilen(sccorname))
2389 write (iout,*) "Bond & inertia constant file : ",
2390 & bondname(:ilen(bondname))
2391 write (iout,*) "Bending parameter file : ",
2392 & thetname(:ilen(thetname))
2393 write (iout,*) "Rotamer parameter file : ",
2394 & rotname(:ilen(rotname))
2395 write (iout,*) "Threading database : ",
2396 & patname(:ilen(patname))
2398 &write (iout,*)" DIRTMP : ",
2400 write (iout,'(80(1h-))')
2404 c----------------------------------------------------------------------------
2405 subroutine card_concat(card)
2406 implicit real*8 (a-h,o-z)
2407 include 'DIMENSIONS'
2408 include 'COMMON.IOUNITS'
2410 character*80 karta,ucase
2412 read (inp,'(a)') karta
2415 do while (karta(80:80).eq.'&')
2416 card=card(:ilen(card)+1)//karta(:79)
2417 read (inp,'(a)') karta
2420 card=card(:ilen(card)+1)//karta
2423 c----------------------------------------------------------------------------------
2425 implicit real*8 (a-h,o-z)
2426 include 'DIMENSIONS'
2427 include 'COMMON.CHAIN'
2428 include 'COMMON.IOUNITS'
2430 include 'COMMON.CONTROL'
2431 open(irest2,file=rest2name,status='unknown')
2432 read(irest2,*) totT,EK,potE,totE,t_bath
2434 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2437 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2439 if(usampl.or.homol_nset.gt.1) then
2440 read (irest2,*) iset
2445 c---------------------------------------------------------------------------------
2446 subroutine read_fragments
2447 implicit real*8 (a-h,o-z)
2448 include 'DIMENSIONS'
2452 include 'COMMON.SETUP'
2453 include 'COMMON.CHAIN'
2454 include 'COMMON.IOUNITS'
2456 include 'COMMON.CONTROL'
2457 read(inp,*) nset,nfrag,npair,nfrag_back
2458 if(me.eq.king.or..not.out1file)
2459 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2460 & " nfrag_back",nfrag_back
2462 read(inp,*) mset(iset)
2464 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2466 if(me.eq.king.or..not.out1file)
2467 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2468 & ifrag(2,i,iset), qinfrag(i,iset)
2471 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2473 if(me.eq.king.or..not.out1file)
2474 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2475 & ipair(2,i,iset), qinpair(i,iset)
2478 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2479 & wfrag_back(3,i,iset),
2480 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2481 if(me.eq.king.or..not.out1file)
2482 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2483 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2488 c-------------------------------------------------------------------------------
2489 subroutine read_dist_constr
2490 implicit real*8 (a-h,o-z)
2491 include 'DIMENSIONS'
2495 include 'COMMON.SETUP'
2496 include 'COMMON.CONTROL'
2497 include 'COMMON.CHAIN'
2498 include 'COMMON.IOUNITS'
2499 include 'COMMON.SBRIDGE'
2500 integer ifrag_(2,100),ipair_(2,100)
2501 double precision wfrag_(100),wpair_(100)
2502 character*500 controlcard
2503 c write (iout,*) "Calling read_dist_constr"
2504 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2506 call card_concat(controlcard)
2507 call readi(controlcard,"NFRAG",nfrag_,0)
2508 call readi(controlcard,"NPAIR",npair_,0)
2509 call readi(controlcard,"NDIST",ndist_,0)
2510 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2511 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2512 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2513 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2514 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2515 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2516 c write (iout,*) "IFRAG"
2518 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2520 c write (iout,*) "IPAIR"
2522 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2524 if (.not.refstr .and. nfrag.gt.0) then
2526 & "ERROR: no reference structure to compute distance restraints"
2528 & "Restraints must be specified explicitly (NDIST=number)"
2531 if (nfrag.lt.2 .and. npair.gt.0) then
2532 write (iout,*) "ERROR: Less than 2 fragments specified",
2533 & " but distance restraints between pairs requested"
2538 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2539 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2540 & ifrag_(2,i)=nstart_sup+nsup-1
2541 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2543 if (wfrag_(i).gt.0.0d0) then
2544 do j=ifrag_(1,i),ifrag_(2,i)-1
2545 do k=j+1,ifrag_(2,i)
2546 c write (iout,*) "j",j," k",k
2548 if (constr_dist.eq.1) then
2553 forcon(nhpb)=wfrag_(i)
2554 else if (constr_dist.eq.2) then
2555 if (ddjk.le.dist_cut) then
2560 forcon(nhpb)=wfrag_(i)
2567 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2570 if (.not.out1file .or. me.eq.king)
2571 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2572 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2574 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2575 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2582 if (wpair_(i).gt.0.0d0) then
2590 do j=ifrag_(1,ii),ifrag_(2,ii)
2591 do k=ifrag_(1,jj),ifrag_(2,jj)
2595 forcon(nhpb)=wpair_(i)
2596 dhpb(nhpb)=dist(j,k)
2598 if (.not.out1file .or. me.eq.king)
2599 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2600 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2602 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2603 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2610 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2611 & ibecarb(i),forcon(nhpb+1)
2612 if (forcon(nhpb+1).gt.0.0d0) then
2614 if (ibecarb(i).gt.0) then
2615 ihpb(i)=ihpb(i)+nres
2616 jhpb(i)=jhpb(i)+nres
2618 if (dhpb(nhpb).eq.0.0d0)
2619 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2623 if (.not.out1file .or. me.eq.king) then
2626 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2627 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2635 c-------------------------------------------------------------------------------
2637 subroutine read_constr_homology
2639 include 'DIMENSIONS'
2643 include 'COMMON.SETUP'
2644 include 'COMMON.CONTROL'
2645 include 'COMMON.CHAIN'
2646 include 'COMMON.IOUNITS'
2648 include 'COMMON.GEO'
2649 include 'COMMON.INTERACT'
2651 c For new homol impl
2653 include 'COMMON.VAR'
2656 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2658 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2659 c & sigma_odl_temp(maxres,maxres,max_template)
2661 character*24 model_ki_dist, model_ki_angle
2662 character*500 controlcard
2663 integer ki, i, j, k, l
2664 logical lprn /.true./
2666 c FP - Nov. 2014 Temporary specifications for new vars
2668 double precision rescore_tmp,x12,y12,z12
2669 double precision, dimension (max_template,maxres) :: rescore
2670 character*24 pdbfile,tpl_k_rescore
2671 c -----------------------------------------------------------------
2672 c Reading multiple PDB ref structures and calculation of retraints
2673 c not using pre-computed ones stored in files model_ki_{dist,angle}
2675 c -----------------------------------------------------------------
2678 c Alternative: reading from input
2679 call card_concat(controlcard)
2680 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2681 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2682 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2683 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2684 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2686 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2687 if (homol_nset.gt.1)then
2688 call card_concat(controlcard)
2689 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2690 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2691 write(iout,*) "iset homology_weight "
2693 write(iout,*) i,waga_homology(i)
2696 iset=mod(kolor,homol_nset)+1
2699 waga_homology(1)=1.0
2702 cd write (iout,*) "nnt",nnt," nct",nct
2714 c Reading HM global scores (prob not required)
2716 c open (4,file="HMscore")
2717 c do k=1,constr_homology
2718 c read (4,*,end=521) hmscore_tmp
2719 c hmscore(k)=hmscore_tmp ! Another transformation can be used
2720 c write(*,*) "Model", k, ":", hmscore(k)
2724 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2725 do k=1,constr_homology
2727 read(inp,'(a)') pdbfile
2728 c Next stament causes error upon compilation (?)
2729 c if(me.eq.king.or. .not. out1file)
2730 c write (iout,'(2a)') 'PDB data will be read from file ',
2731 c & pdbfile(:ilen(pdbfile))
2732 open(ipdbin,file=pdbfile,status='old',err=33)
2734 33 write (iout,'(a)') 'Error opening PDB file.'
2737 c print *,'Begin reading pdb data'
2739 c Files containing res sim or local scores (former containing sigmas)
2742 write(kic2,'(bz,i2.2)') k
2744 tpl_k_rescore="template"//kic2//".sco"
2745 c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
2746 c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
2747 c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
2748 c tpl_k_sigma_d="template"//kic2//".sigma_d"
2753 c Distance restraints
2756 C Copy the coordinates from reference coordinates (?)
2760 c write (iout,*) "c(",j,i,") =",c(j,i)
2764 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2766 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2767 open (ientin,file=tpl_k_rescore,status='old')
2768 do irec=1,maxdim ! loop for reading res sim
2770 rescore(k,irec)=0.0d0
2773 read (ientin,*,end=1401) rescore_tmp
2774 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2775 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2776 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2781 c open (ientin,file=tpl_k_sigma_odl,status='old')
2782 c do irec=1,maxdim ! loop for reading sigma_odl
2783 c read (ientin,*,end=1401) i, j,
2784 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
2785 c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
2786 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
2790 if (waga_dist.gt.0.0d0) then
2792 do i = nnt,nct-2 ! right? without parallel.
2793 do j=i+2,nct ! right?
2794 c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
2795 c do j=i+2,nres ! ibid
2796 c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
2797 c do j=i+2,nct ! ibid
2799 c write (iout,*) "k",k
2800 c write (iout,*) "i",i," j",j," constr_homology",
2805 c Attempt to replace dist(i,j) by its definition in ...
2810 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2813 c odl(k,ii)=dist(i,j)
2814 c write (iout,*) "dist(",i,j,") =",dist(i,j)
2815 c write (iout,*) "distal = ",distal
2816 c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
2817 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2818 c & "rescore(",k,j,") =",rescore(k,j)
2820 c Calculation of sigma from res sim
2822 c if (odl(k,ii).le.6.0d0) then
2823 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2824 c Other functional forms possible depending on odl(k,ii), eg.
2826 if (odl(k,ii).le.dist_cut) then
2827 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
2828 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2830 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
2831 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2833 c Following expr replaced by a positive exp argument
2834 c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2835 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2837 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
2838 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2841 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
2842 c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
2844 c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
2845 c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
2847 c read (ientin,*) sigma_odl(k,ii) ! 1st variant
2850 c if (constr_homology.gt.0) call homology_partition
2853 c Theta, dihedral and SC retraints
2855 if (waga_angle.gt.0.0d0) then
2856 c open (ientin,file=tpl_k_sigma_dih,status='old')
2857 c do irec=1,maxres-3 ! loop for reading sigma_dih
2858 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2859 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2860 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2861 c & sigma_dih(k,i+nnt-1)
2865 do i = nnt+3,nct ! right? without parallel.
2866 c do i=1,nres ! alternative for bounds acc to readpdb?
2867 c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
2868 c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
2869 dih(k,i)=phiref(i) ! right?
2870 c read (ientin,*) sigma_dih(k,i) ! original variant
2871 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2872 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2873 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2874 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2875 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2877 sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
2878 & rescore(k,i-2)+rescore(k,i-3) ! right expression ?
2880 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2881 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2882 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2883 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2884 c Instead of res sim other local measure of b/b str reliability possible
2885 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2886 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2887 if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
2888 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
2892 if (waga_theta.gt.0.0d0) then
2893 c open (ientin,file=tpl_k_sigma_theta,status='old')
2894 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2895 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2896 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2897 c & sigma_theta(k,i+nnt-1)
2902 do i = nnt+2,nct ! right? without parallel.
2903 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2904 c do i=ithet_start,ithet_end ! with FG parallel.
2905 thetatpl(k,i)=thetaref(i)
2906 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2907 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2908 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2909 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2910 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2911 sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
2912 & rescore(k,i-2) ! right expression ?
2913 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2915 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2916 c rescore(k,i-2) ! right expression ?
2917 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2918 if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
2922 if (waga_d.gt.0.0d0) then
2923 c open (ientin,file=tpl_k_sigma_d,status='old')
2924 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2925 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2926 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2927 c & sigma_d(k,i+nnt-1)
2932 do i = nnt,nct ! right? without parallel.
2933 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2934 c do i=loc_start,loc_end ! with FG parallel.
2935 if (itype(i).eq.10) goto 1 ! right?
2939 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2940 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2941 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2942 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
2943 sigma_d(k,i)=rescore(k,i) ! right expression ?
2944 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2946 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
2947 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2948 c read (ientin,*) sigma_d(k,i) ! 1st variant
2949 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
2955 if (waga_dist.gt.0.0d0) lim_odl=ii
2956 if (constr_homology.gt.0) call homology_partition
2957 if (constr_homology.gt.0) call init_int_table
2958 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
2959 cd & "lim_xx=",lim_xx
2960 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2961 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2965 if (.not.lprn) return
2966 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2967 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2968 write (iout,*) "Distance restraints from templates"
2970 write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
2971 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
2973 write (iout,*) "Dihedral angle restraints from templates"
2975 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
2976 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2978 write (iout,*) "Virtual-bond angle restraints from templates"
2979 do i=nnt+2,lim_theta
2980 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
2981 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2983 write (iout,*) "SC restraints from templates"
2985 write(iout,'(i5,10(4f8.2,4x))') i,
2986 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
2987 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2990 c -----------------------------------------------------------------
2993 c----------------------------------------------------------------------
2996 subroutine flush(iu)
3001 subroutine flush(iu)
3007 c------------------------------------------------------------------------------
3008 subroutine copy_to_tmp(source)
3009 include "DIMENSIONS"
3010 include "COMMON.IOUNITS"
3011 character*(*) source
3012 character* 256 tmpfile
3016 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3017 inquire(file=tmpfile,exist=ex)
3019 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3020 & " to temporary directory..."
3021 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3022 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3026 c------------------------------------------------------------------------------
3027 subroutine move_from_tmp(source)
3028 include "DIMENSIONS"
3029 include "COMMON.IOUNITS"
3030 character*(*) source
3033 write (*,*) "Moving ",source(:ilen(source)),
3034 & " from temporary directory to working directory"
3035 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3036 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3039 c------------------------------------------------------------------------------
3040 subroutine random_init(seed)
3042 C Initialize random number generator
3044 implicit real*8 (a-h,o-z)
3045 include 'DIMENSIONS'
3051 logical OKRandom, prng_restart
3053 integer iseed_array(4)
3055 include 'COMMON.IOUNITS'
3056 include 'COMMON.TIME1'
3057 include 'COMMON.THREAD'
3058 include 'COMMON.SBRIDGE'
3059 include 'COMMON.CONTROL'
3060 include 'COMMON.MCM'
3061 include 'COMMON.MAP'
3062 include 'COMMON.HEADER'
3063 csa include 'COMMON.CSA'
3064 include 'COMMON.CHAIN'
3065 include 'COMMON.MUCA'
3067 include 'COMMON.FFIELD'
3068 include 'COMMON.SETUP'
3069 iseed=-dint(dabs(seed))
3070 if (iseed.eq.0) then
3071 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3072 & 'Random seed undefined. The program will stop.'
3073 write (*,'(/80(1h*)/20x,a/80(1h*))')
3074 & 'Random seed undefined. The program will stop.'
3076 call mpi_finalize(mpi_comm_world,ierr)
3078 stop 'Bad random seed.'
3081 if (fg_rank.eq.0) then
3085 if(me.eq.king .or. .not. out1file)
3086 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3087 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3088 OKRandom = prng_restart(me,iseedi8)
3091 tmp=65536.0d0**(4-i)
3092 iseed_array(i) = dint(seed/tmp)
3093 seed=seed-iseed_array(i)*tmp
3095 if(me.eq.king .or. .not. out1file)
3096 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3097 & (iseed_array(i),i=1,4)
3098 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3099 & (iseed_array(i),i=1,4)
3100 OKRandom = prng_restart(me,iseed_array)
3103 r1=ran_number(0.0D0,1.0D0)
3104 if(me.eq.king .or. .not. out1file)
3105 & write (iout,*) 'ran_num',r1
3106 if (r1.lt.0.0d0) OKRandom=.false.
3108 if (.not.OKRandom) then
3109 write (iout,*) 'PRNG IS NOT WORKING!!!'
3110 print *,'PRNG IS NOT WORKING!!!'
3113 call mpi_abort(mpi_comm_world,error_msg,ierr)
3116 write (iout,*) 'too many processors for parallel prng'
3117 write (*,*) 'too many processors for parallel prng'
3125 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)