2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(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'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
102 call readi(controlcard,'SYM',symetr,1)
103 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
104 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
105 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
106 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
107 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
108 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
109 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
110 call reada(controlcard,'DRMS',drms,0.1D0)
111 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
112 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
113 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
114 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
115 write (iout,'(a,f10.1)')'DRMS = ',drms
116 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
117 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
119 call readi(controlcard,'NZ_START',nz_start,0)
120 call readi(controlcard,'NZ_END',nz_end,0)
121 call readi(controlcard,'IZ_SC',iz_sc,0)
123 safety = 60.0d0*safety
126 call reada(controlcard,"T_BATH",t_bath,300.0d0)
127 minim=(index(controlcard,'MINIMIZE').gt.0)
128 dccart=(index(controlcard,'CART').gt.0)
129 overlapsc=(index(controlcard,'OVERLAP').gt.0)
130 overlapsc=.not.overlapsc
131 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
132 searchsc=.not.searchsc
133 sideadd=(index(controlcard,'SIDEADD').gt.0)
134 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
135 outpdb=(index(controlcard,'PDBOUT').gt.0)
136 outmol2=(index(controlcard,'MOL2OUT').gt.0)
137 pdbref=(index(controlcard,'PDBREF').gt.0)
138 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
139 indpdb=index(controlcard,'PDBSTART')
140 extconf=(index(controlcard,'EXTCONF').gt.0)
141 AFMlog=(index(controlcard,'AFM'))
142 selfguide=(index(controlcard,'SELFGUIDE'))
143 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
144 call readi(controlcard,'IPRINT',iprint,0)
145 call readi(controlcard,'MAXGEN',maxgen,10000)
146 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
147 call readi(controlcard,"KDIAG",kdiag,0)
148 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
149 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
150 & write (iout,*) "RESCALE_MODE",rescale_mode
151 split_ene=index(controlcard,'SPLIT_ENE').gt.0
152 if (index(controlcard,'REGULAR').gt.0.0D0) then
153 call reada(controlcard,'WEIDIS',weidis,0.1D0)
157 if (index(controlcard,'CHECKGRAD').gt.0) then
159 if (index(controlcard,'CART').gt.0) then
161 elseif (index(controlcard,'CARINT').gt.0) then
166 elseif (index(controlcard,'THREAD').gt.0) then
168 call readi(controlcard,'THREAD',nthread,0)
169 if (nthread.gt.0) then
170 call reada(controlcard,'WEIDIS',weidis,0.1D0)
173 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
174 stop 'Error termination in Read_Control.'
176 else if (index(controlcard,'MCMA').gt.0) then
178 else if (index(controlcard,'MCEE').gt.0) then
180 else if (index(controlcard,'MULTCONF').gt.0) then
182 else if (index(controlcard,'MAP').gt.0) then
184 call readi(controlcard,'MAP',nmap,0)
185 else if (index(controlcard,'CSA').gt.0) then
187 crc else if (index(controlcard,'ZSCORE').gt.0) then
189 crc ZSCORE is rm from UNRES, modecalc=9 is available
192 cfcm else if (index(controlcard,'MCMF').gt.0) then
194 else if (index(controlcard,'SOFTREG').gt.0) then
196 else if (index(controlcard,'CHECK_BOND').gt.0) then
198 else if (index(controlcard,'TEST').gt.0) then
200 else if (index(controlcard,'MD').gt.0) then
202 else if (index(controlcard,'RE ').gt.0) then
206 lmuca=index(controlcard,'MUCA').gt.0
207 call readi(controlcard,'MUCADYN',mucadyn,0)
208 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
209 if (lmuca .and. (me.eq.king .or. .not.out1file ))
211 write (iout,*) 'MUCADYN=',mucadyn
212 write (iout,*) 'MUCASMOOTH=',muca_smooth
215 iscode=index(controlcard,'ONE_LETTER')
216 indphi=index(controlcard,'PHI')
217 indback=index(controlcard,'BACK')
218 iranconf=index(controlcard,'RAND_CONF')
219 i2ndstr=index(controlcard,'USE_SEC_PRED')
220 gradout=index(controlcard,'GRADOUT').gt.0
221 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
222 C DISTCHAINMAX become obsolete for periodic boundry condition
223 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
224 C Reading the dimensions of box in x,y,z coordinates
225 call reada(controlcard,'BOXX',boxxsize,100.0d0)
226 call reada(controlcard,'BOXY',boxysize,100.0d0)
227 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
228 c Cutoff range for interactions
229 call reada(controlcard,"R_CUT",r_cut,15.0d0)
230 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
231 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
232 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
233 if (lipthick.gt.0.0d0) then
234 bordliptop=(boxzsize+lipthick)/2.0
235 bordlipbot=bordliptop-lipthick
237 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
238 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
239 buflipbot=bordlipbot+lipbufthick
240 bufliptop=bordliptop-lipbufthick
241 if ((lipbufthick*2.0d0).gt.lipthick)
242 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
244 write(iout,*) "bordliptop=",bordliptop
245 write(iout,*) "bordlipbot=",bordlipbot
246 write(iout,*) "bufliptop=",bufliptop
247 write(iout,*) "buflipbot=",buflipbot
250 if (me.eq.king .or. .not.out1file )
251 & write (iout,*) "DISTCHAINMAX",distchainmax
253 if(me.eq.king.or..not.out1file)
254 & write (iout,'(2a)') diagmeth(kdiag),
255 & ' routine used to diagonalize matrices.'
258 c--------------------------------------------------------------------------
259 subroutine read_REMDpar
263 implicit real*8 (a-h,o-z)
265 include 'COMMON.IOUNITS'
266 include 'COMMON.TIME1'
269 include 'COMMON.LANGEVIN'
271 include 'COMMON.LANGEVIN.lang0'
273 include 'COMMON.INTERACT'
274 include 'COMMON.NAMES'
276 include 'COMMON.REMD'
277 include 'COMMON.CONTROL'
278 include 'COMMON.SETUP'
280 character*320 controlcard
281 character*3200 controlcard1
282 integer iremd_m_total
284 if(me.eq.king.or..not.out1file)
285 & write (iout,*) "REMD setup"
287 call card_concat(controlcard)
288 call readi(controlcard,"NREP",nrep,3)
289 call readi(controlcard,"NSTEX",nstex,1000)
290 call reada(controlcard,"RETMIN",retmin,10.0d0)
291 call reada(controlcard,"RETMAX",retmax,1000.0d0)
292 mremdsync=(index(controlcard,'SYNC').gt.0)
293 call readi(controlcard,"NSYN",i_sync_step,100)
294 restart1file=(index(controlcard,'REST1FILE').gt.0)
295 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
296 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
297 if(max_cache_traj_use.gt.max_cache_traj)
298 & max_cache_traj_use=max_cache_traj
299 if(me.eq.king.or..not.out1file) then
300 cd if (traj1file) then
301 crc caching is in testing - NTWX is not ignored
302 cd write (iout,*) "NTWX value is ignored"
303 cd write (iout,*) " trajectory is stored to one file by master"
304 cd write (iout,*) " before exchange at NSTEX intervals"
306 write (iout,*) "NREP= ",nrep
307 write (iout,*) "NSTEX= ",nstex
308 write (iout,*) "SYNC= ",mremdsync
309 write (iout,*) "NSYN= ",i_sync_step
310 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
313 if (index(controlcard,'TLIST').gt.0) then
315 call card_concat(controlcard1)
316 read(controlcard1,*) (remd_t(i),i=1,nrep)
317 if(me.eq.king.or..not.out1file)
318 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
321 if (index(controlcard,'MLIST').gt.0) then
323 call card_concat(controlcard1)
324 read(controlcard1,*) (remd_m(i),i=1,nrep)
325 if(me.eq.king.or..not.out1file) then
326 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
329 iremd_m_total=iremd_m_total+remd_m(i)
331 write (iout,*) 'Total number of replicas ',iremd_m_total
334 if(me.eq.king.or..not.out1file)
335 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
338 c--------------------------------------------------------------------------
339 subroutine read_MDpar
343 implicit real*8 (a-h,o-z)
345 include 'COMMON.IOUNITS'
346 include 'COMMON.TIME1'
349 include 'COMMON.LANGEVIN'
351 include 'COMMON.LANGEVIN.lang0'
353 include 'COMMON.INTERACT'
354 include 'COMMON.NAMES'
356 include 'COMMON.SETUP'
357 include 'COMMON.CONTROL'
358 include 'COMMON.SPLITELE'
360 character*320 controlcard
362 call card_concat(controlcard)
363 call readi(controlcard,"NSTEP",n_timestep,1000000)
364 call readi(controlcard,"NTWE",ntwe,100)
365 call readi(controlcard,"NTWX",ntwx,1000)
366 call reada(controlcard,"DT",d_time,1.0d-1)
367 call reada(controlcard,"DVMAX",dvmax,2.0d1)
368 call reada(controlcard,"DAMAX",damax,1.0d1)
369 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
370 call readi(controlcard,"LANG",lang,0)
371 RESPA = index(controlcard,"RESPA") .gt. 0
372 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
373 ntime_split0=ntime_split
374 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
375 ntime_split0=ntime_split
376 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
377 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
378 rest = index(controlcard,"REST").gt.0
379 tbf = index(controlcard,"TBF").gt.0
380 usampl = index(controlcard,"USAMPL").gt.0
382 mdpdb = index(controlcard,"MDPDB").gt.0
383 call reada(controlcard,"T_BATH",t_bath,300.0d0)
384 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
385 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
386 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
387 if (count_reset_moment.eq.0) count_reset_moment=1000000000
388 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
389 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
390 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
391 if (count_reset_vel.eq.0) count_reset_vel=1000000000
392 large = index(controlcard,"LARGE").gt.0
393 print_compon = index(controlcard,"PRINT_COMPON").gt.0
394 rattle = index(controlcard,"RATTLE").gt.0
395 preminim = index(controlcard,"PREMINIM").gt.0
397 dccart=(index(controlcard,'CART').gt.0)
400 c if performing umbrella sampling, fragments constrained are read from the fragment file
406 if(me.eq.king.or..not.out1file) then
408 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
410 write (iout,'(a)') "The units are:"
411 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
412 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
413 & " acceleration: angstrom/(48.9 fs)**2"
414 write (iout,'(a)') "energy: kcal/mol, temperature: K"
416 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
417 write (iout,'(a60,f10.5,a)')
418 & "Initial time step of numerical integration:",d_time,
420 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
422 write (iout,'(2a,i4,a)')
423 & "A-MTS algorithm used; initial time step for fast-varying",
424 & " short-range forces split into",ntime_split," steps."
425 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
426 & r_cut," lambda",rlamb
428 write (iout,'(2a,f10.5)')
429 & "Maximum acceleration threshold to reduce the time step",
430 & "/increase split number:",damax
431 write (iout,'(2a,f10.5)')
432 & "Maximum predicted energy drift to reduce the timestep",
433 & "/increase split number:",edriftmax
434 write (iout,'(a60,f10.5)')
435 & "Maximum velocity threshold to reduce velocities:",dvmax
436 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
437 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
438 if (rattle) write (iout,'(a60)')
439 & "Rattle algorithm used to constrain the virtual bonds"
440 if (preminim .or. iranconf.gt.0) then
442 & "Initial structure will be energy-minimized"
447 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
448 call reada(controlcard,"RWAT",rwat,1.4d0)
449 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
450 surfarea=index(controlcard,"SURFAREA").gt.0
451 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
452 if(me.eq.king.or..not.out1file)then
453 write (iout,'(/a,$)') "Langevin dynamics calculation"
456 & " with direct integration of Langevin equations"
457 else if (lang.eq.2) then
458 write (iout,'(a/)') " with TINKER stochasic MD integrator"
459 else if (lang.eq.3) then
460 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
461 else if (lang.eq.4) then
462 write (iout,'(a/)') " in overdamped mode"
464 write (iout,'(//a,i5)')
465 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
468 write (iout,'(a60,f10.5)') "Temperature:",t_bath
469 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
470 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
471 write (iout,'(a60,f10.5)')
472 & "Scaling factor of the friction forces:",scal_fric
473 if (surfarea) write (iout,'(2a,i10,a)')
474 & "Friction coefficients will be scaled by solvent-accessible",
475 & " surface area every",reset_fricmat," steps."
477 c Calculate friction coefficients and bounds of stochastic forces
478 eta=6*pi*cPoise*etawat
479 if(me.eq.king.or..not.out1file)
480 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
482 gamp=scal_fric*(pstok+rwat)*eta
483 stdfp=dsqrt(2*Rb*t_bath/d_time)
485 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
486 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
488 if(me.eq.king.or..not.out1file)then
489 write (iout,'(/2a/)')
490 & "Radii of site types and friction coefficients and std's of",
491 & " stochastic forces of fully exposed sites"
492 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
494 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
495 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
499 if(me.eq.king.or..not.out1file)then
500 write (iout,'(a)') "Berendsen bath calculation"
501 write (iout,'(a60,f10.5)') "Temperature:",t_bath
502 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
504 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
505 & count_reset_moment," steps"
507 & write (iout,'(a,i10,a)')
508 & "Velocities will be reset at random every",count_reset_vel,
512 if(me.eq.king.or..not.out1file)
513 & write (iout,'(a31)') "Microcanonical mode calculation"
515 if(me.eq.king.or..not.out1file)then
516 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
518 write(iout,*) "MD running with constraints."
519 write(iout,*) "Equilibration time ", eq_time, " mtus."
520 write(iout,*) "Constraining ", nfrag," fragments."
521 write(iout,*) "Length of each fragment, weight and q0:"
523 write (iout,*) "Set of restraints #",iset
525 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
526 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
528 write(iout,*) "constraints between ", npair, "fragments."
529 write(iout,*) "constraint pairs, weights and q0:"
531 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
532 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
534 write(iout,*) "angle constraints within ", nfrag_back,
535 & "backbone fragments."
536 write(iout,*) "fragment, weights:"
538 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
539 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
540 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
543 iset=mod(kolor,nset)+1
546 if(me.eq.king.or..not.out1file)
547 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
550 c------------------------------------------------------------------------------
553 C Read molecular data.
555 implicit real*8 (a-h,o-z)
561 include 'COMMON.IOUNITS'
564 include 'COMMON.INTERACT'
565 include 'COMMON.LOCAL'
566 include 'COMMON.NAMES'
567 include 'COMMON.CHAIN'
568 include 'COMMON.FFIELD'
569 include 'COMMON.SBRIDGE'
570 include 'COMMON.HEADER'
571 include 'COMMON.CONTROL'
572 include 'COMMON.DBASE'
573 include 'COMMON.THREAD'
574 include 'COMMON.CONTACTS'
575 include 'COMMON.TORCNSTR'
576 include 'COMMON.TIME1'
577 include 'COMMON.BOUNDS'
579 include 'COMMON.SETUP'
580 character*4 sequence(maxres)
582 double precision x(maxvar)
583 character*256 pdbfile
584 character*320 weightcard
585 character*80 weightcard_t,ucase
586 dimension itype_pdb(maxres)
587 common /pizda/ itype_pdb
588 logical seq_comp,fail
589 double precision energia(0:n_ene)
595 C Read weights of the subsequent energy terms.
596 call card_concat(weightcard)
597 call reada(weightcard,'WLONG',wlong,1.0D0)
598 call reada(weightcard,'WSC',wsc,wlong)
599 call reada(weightcard,'WSCP',wscp,wlong)
600 call reada(weightcard,'WELEC',welec,1.0D0)
601 call reada(weightcard,'WVDWPP',wvdwpp,welec)
602 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
603 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
604 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
605 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
606 call reada(weightcard,'WTURN3',wturn3,1.0D0)
607 call reada(weightcard,'WTURN4',wturn4,1.0D0)
608 call reada(weightcard,'WTURN6',wturn6,1.0D0)
609 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
610 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
611 call reada(weightcard,'WBOND',wbond,1.0D0)
612 call reada(weightcard,'WTOR',wtor,1.0D0)
613 call reada(weightcard,'WTORD',wtor_d,1.0D0)
614 call reada(weightcard,'WANG',wang,1.0D0)
615 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
616 call reada(weightcard,'SCAL14',scal14,0.4D0)
617 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
618 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
619 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
620 call reada(weightcard,'TEMP0',temp0,300.0d0)
621 call reada(weightcard,'WLT',wliptran,0.0D0)
622 if (index(weightcard,'SOFT').gt.0) ipot=6
623 C 12/1/95 Added weight for the multi-body term WCORR
624 call reada(weightcard,'WCORRH',wcorr,1.0D0)
625 if (wcorr4.gt.0.0d0) wcorr=wcorr4
645 if(me.eq.king.or..not.out1file)
646 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
647 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
649 10 format (/'Energy-term weights (unscaled):'//
650 & 'WSCC= ',f10.6,' (SC-SC)'/
651 & 'WSCP= ',f10.6,' (SC-p)'/
652 & 'WELEC= ',f10.6,' (p-p electr)'/
653 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
654 & 'WBOND= ',f10.6,' (stretching)'/
655 & 'WANG= ',f10.6,' (bending)'/
656 & 'WSCLOC= ',f10.6,' (SC local)'/
657 & 'WTOR= ',f10.6,' (torsional)'/
658 & 'WTORD= ',f10.6,' (double torsional)'/
659 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
660 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
661 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
662 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
663 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
664 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
665 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
666 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
667 & 'WTURN6= ',f10.6,' (turns, 6th order)')
668 if(me.eq.king.or..not.out1file)then
669 if (wcorr4.gt.0.0d0) then
670 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
671 & 'between contact pairs of peptide groups'
672 write (iout,'(2(a,f5.3/))')
673 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
674 & 'Range of quenching the correlation terms:',2*delt_corr
675 else if (wcorr.gt.0.0d0) then
676 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
677 & 'between contact pairs of peptide groups'
679 write (iout,'(a,f8.3)')
680 & 'Scaling factor of 1,4 SC-p interactions:',scal14
681 write (iout,'(a,f8.3)')
682 & 'General scaling factor of SC-p interactions:',scalscp
684 r0_corr=cutoff_corr-delt_corr
686 aad(i,1)=scalscp*aad(i,1)
687 aad(i,2)=scalscp*aad(i,2)
688 bad(i,1)=scalscp*bad(i,1)
689 bad(i,2)=scalscp*bad(i,2)
691 call rescale_weights(t_bath)
692 if(me.eq.king.or..not.out1file)
693 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
694 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
696 22 format (/'Energy-term weights (scaled):'//
697 & 'WSCC= ',f10.6,' (SC-SC)'/
698 & 'WSCP= ',f10.6,' (SC-p)'/
699 & 'WELEC= ',f10.6,' (p-p electr)'/
700 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
701 & 'WBOND= ',f10.6,' (stretching)'/
702 & 'WANG= ',f10.6,' (bending)'/
703 & 'WSCLOC= ',f10.6,' (SC local)'/
704 & 'WTOR= ',f10.6,' (torsional)'/
705 & 'WTORD= ',f10.6,' (double torsional)'/
706 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
707 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
708 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
709 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
710 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
711 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
712 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
713 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
714 & 'WTURN6= ',f10.6,' (turns, 6th order)')
715 if(me.eq.king.or..not.out1file)
716 & write (iout,*) "Reference temperature for weights calculation:",
718 call reada(weightcard,"D0CM",d0cm,3.78d0)
719 call reada(weightcard,"AKCM",akcm,15.1d0)
720 call reada(weightcard,"AKTH",akth,11.0d0)
721 call reada(weightcard,"AKCT",akct,12.0d0)
722 call reada(weightcard,"V1SS",v1ss,-1.08d0)
723 call reada(weightcard,"V2SS",v2ss,7.61d0)
724 call reada(weightcard,"V3SS",v3ss,13.7d0)
725 call reada(weightcard,"EBR",ebr,-5.50D0)
726 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
728 dyn_ss_mask(i)=.false.
732 dyn_ssbond_ij(i,j)=1.0d300
735 call reada(weightcard,"HT",Ht,0.0D0)
737 ss_depth=ebr/wsc-0.25*eps(1,1)
738 Ht=Ht/wsc-0.25*eps(1,1)
739 akcm=akcm*wstrain/wsc
740 akth=akth*wstrain/wsc
741 akct=akct*wstrain/wsc
742 v1ss=v1ss*wstrain/wsc
743 v2ss=v2ss*wstrain/wsc
744 v3ss=v3ss*wstrain/wsc
746 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
749 if(me.eq.king.or..not.out1file) then
750 write (iout,*) "Parameters of the SS-bond potential:"
751 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
753 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
754 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
755 write (iout,*)" HT",Ht
756 print *,'indpdb=',indpdb,' pdbref=',pdbref
758 if (indpdb.gt.0 .or. pdbref) then
759 read(inp,'(a)') pdbfile
760 if(me.eq.king.or..not.out1file)
761 & write (iout,'(2a)') 'PDB data will be read from file ',
762 & pdbfile(:ilen(pdbfile))
763 open(ipdbin,file=pdbfile,status='old',err=33)
765 33 write (iout,'(a)') 'Error opening PDB file.'
768 c print *,'Begin reading pdb data'
777 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
778 & (crefjlee(j,i+nres),j=1,3)
781 c print *,'Finished reading pdb data'
782 if(me.eq.king.or..not.out1file)
783 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
784 & ' nstart_sup=',nstart_sup
786 itype_pdb(i)=itype(i)
790 nct=nstart_sup+nsup-1
791 call contact(.false.,ncont_ref,icont_ref,co)
794 if(me.eq.king.or..not.out1file)
795 & write(iout,*)'Adding sidechains'
799 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
802 do while (fail.and.nsi.le.maxsi)
803 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
806 if(fail) write(iout,*)'Adding sidechain failed for res ',
807 & i,' after ',nsi,' trials'
810 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
814 if (indpdb.eq.0) then
815 C Read sequence if not taken from the pdb file.
817 c print *,'nres=',nres
818 if (iscode.gt.0) then
819 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
821 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
823 C Convert sequence to numeric code
825 itype(i)=rescode(i,sequence(i),iscode)
827 C Assign initial virtual bond lengths
833 vbld(i+nres)=dsc(iabs(itype(i)))
834 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
835 c write (iout,*) "i",i," itype",itype(i),
836 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
840 c print '(20i4)',(itype(i),i=1,nres)
843 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
845 if (itype(i).eq.ntyp1) then
849 else if (iabs(itype(i+1)).ne.20) then
851 else if (iabs(itype(i)).ne.20) then
858 if(me.eq.king.or..not.out1file)then
859 write (iout,*) "ITEL"
861 write (iout,*) i,itype(i),itel(i)
863 print *,'Call Read_Bridge.'
866 C 8/13/98 Set limits to generating the dihedral angles
871 read (inp,*) ndih_constr
872 if (ndih_constr.gt.0) then
874 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
875 if(me.eq.king.or..not.out1file)then
877 & 'There are',ndih_constr,' constraints on phi angles.'
879 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
883 phi0(i)=deg2rad*phi0(i)
884 drange(i)=deg2rad*drange(i)
886 if(me.eq.king.or..not.out1file)
887 & write (iout,*) 'FTORS',ftors
890 phibound(1,ii) = phi0(i)-drange(i)
891 phibound(2,ii) = phi0(i)+drange(i)
898 write (iout,'(a)') 'Boundaries in phi angle sampling:'
900 write (iout,'(a3,i5,2f10.1)')
901 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
907 cd print *,'NNT=',NNT,' NCT=',NCT
908 if (itype(1).eq.ntyp1) nnt=2
909 if (itype(nres).eq.ntyp1) nct=nct-1
911 if(me.eq.king.or..not.out1file)
912 & write (iout,'(a,i3)') 'nsup=',nsup
914 if (nsup.le.(nct-nnt+1)) then
915 do i=0,nct-nnt+1-nsup
916 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
922 & 'Error - sequences to be superposed do not match.'
925 do i=0,nsup-(nct-nnt+1)
926 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
928 nstart_sup=nstart_sup+i
934 & 'Error - sequences to be superposed do not match.'
937 if (nsup.eq.0) nsup=nct-nnt
938 if (nstart_sup.eq.0) nstart_sup=nnt
939 if (nstart_seq.eq.0) nstart_seq=nnt
940 if(me.eq.king.or..not.out1file)
941 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
942 & ' nstart_seq=',nstart_seq
944 c--- Zscore rms -------
945 if (nz_start.eq.0) nz_start=nnt
946 if (nz_end.eq.0 .and. nsup.gt.0) then
948 else if (nz_end.eq.0) then
951 if(me.eq.king.or..not.out1file)then
952 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
953 write (iout,*) 'IZ_SC=',iz_sc
955 c----------------------
958 if (.not.pdbref) then
959 call read_angles(inp,*38)
961 38 write (iout,'(a)') 'Error reading reference structure.'
963 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
964 stop 'Error reading reference structure'
968 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
978 call contact(.true.,ncont_ref,icont_ref,co)
980 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
982 if (constr_dist.gt.0) call read_dist_constr
983 write (iout,*) "After read_dist_constr nhpb",nhpb
984 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
985 if(me.eq.king.or..not.out1file)
986 & write (iout,*) 'Contact order:',co
988 if(me.eq.king.or..not.out1file)
989 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
992 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
994 if(me.eq.king.or..not.out1file)
995 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
996 & icont_ref(1,i),' ',
997 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1003 if (constr_homology.gt.0) then
1004 call read_constr_homology
1005 if (indpdb.gt.0 .or. pdbref) then
1008 c(j,i)=crefjlee(j,i)
1009 cref(j,i,1)=crefjlee(j,i)
1014 write (iout,*) "Array C"
1016 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1017 & (c(j,i+nres),j=1,3)
1019 write (iout,*) "Array Cref"
1021 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1022 & (cref(j,i+nres,1),j=1,3)
1025 call int_from_cart1(.false.)
1026 call sc_loc_geom(.false.)
1028 thetaref(i)=theta(i)
1033 dc(j,i)=c(j,i+1)-c(j,i)
1034 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1039 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1040 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1048 if (nhpb.gt.0) call hpb_partition
1049 c write (iout,*) "After read_dist_constr nhpb",nhpb
1051 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1052 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1053 & modecalc.ne.10) then
1054 C If input structure hasn't been supplied from the PDB file read or generate
1056 if (iranconf.eq.0 .and. .not. extconf) then
1057 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1058 & write (iout,'(a)') 'Initial geometry will be read in.'
1060 read(inp,'(8f10.5)',end=36,err=36)
1061 & ((c(l,k),l=1,3),k=1,nres),
1062 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1063 write (iout,*) "Exit READ_CART"
1064 write (iout,'(8f10.5)')
1065 & ((c(l,k),l=1,3),k=1,nres),
1066 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1067 call int_from_cart1(.true.)
1068 write (iout,*) "Finish INT_TO_CART"
1071 dc(j,i)=c(j,i+1)-c(j,i)
1072 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1076 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1078 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1079 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1085 call read_angles(inp,*36)
1088 36 write (iout,'(a)') 'Error reading angle file.'
1090 call mpi_finalize( MPI_COMM_WORLD,IERR )
1092 stop 'Error reading angle file.'
1094 else if (extconf) then
1095 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1096 & write (iout,'(a)') 'Extended chain initial geometry.'
1098 theta(i)=90d0*deg2rad
1101 phi(i)=180d0*deg2rad
1104 alph(i)=110d0*deg2rad
1107 omeg(i)=-120d0*deg2rad
1108 if (itype(i).le.0) omeg(i)=-omeg(i)
1111 if(me.eq.king.or..not.out1file)
1112 & write (iout,'(a)') 'Random-generated initial geometry.'
1116 if (me.eq.king .or. fg_rank.eq.0 .and. (
1117 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1121 call gen_rand_conf(itmp,*30)
1123 30 write (iout,*) 'Failed to generate random conformation',
1124 & ', itrial=',itrial
1125 write (*,*) 'Processor:',me,
1126 & ' Failed to generate random conformation',
1136 write (iout,'(a,i3,a)') 'Processor:',me,
1137 & ' error in generating random conformation.'
1138 write (*,'(a,i3,a)') 'Processor:',me,
1139 & ' error in generating random conformation.'
1142 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1147 & ' error in generating random conformation.'
1152 elseif (modecalc.eq.4) then
1153 read (inp,'(a)') intinname
1154 open (intin,file=intinname,status='old',err=333)
1155 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1156 & write (iout,'(a)') 'intinname',intinname
1157 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1159 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1161 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1163 stop 'Error opening angle file.'
1167 C Generate distance constraints, if the PDB structure is to be regularized.
1168 if (nthread.gt.0) then
1169 call read_threadbase
1172 if (me.eq.king .or. .not. out1file)
1174 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1175 write (iout,'(/a,i3,a)')
1176 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1177 write (iout,'(20i4)') (iss(i),i=1,ns)
1179 write(iout,*)"Running with dynamic disulfide-bond formation"
1181 write (iout,'(/a/)') 'Pre-formed links are:'
1187 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1188 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1194 if (ns.gt.0.and.dyn_ss) then
1198 forcon(i-nss)=forcon(i)
1205 dyn_ss_mask(iss(i))=.true.
1208 if (i2ndstr.gt.0) call secstrp2dihc
1209 c call geom_to_var(nvar,x)
1210 c call etotal(energia(0))
1211 c call enerprint(energia(0))
1212 c call briefout(0,etot)
1214 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1215 cd write (iout,'(a)') 'Variable list:'
1216 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1218 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1219 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1220 & 'Processor',myrank,': end reading molecular data.'
1224 c--------------------------------------------------------------------------
1225 logical function seq_comp(itypea,itypeb,length)
1227 integer length,itypea(length),itypeb(length)
1230 if (itypea(i).ne.itypeb(i)) then
1238 c-----------------------------------------------------------------------------
1239 subroutine read_bridge
1240 C Read information about disulfide bridges.
1241 implicit real*8 (a-h,o-z)
1242 include 'DIMENSIONS'
1246 include 'COMMON.IOUNITS'
1247 include 'COMMON.GEO'
1248 include 'COMMON.VAR'
1249 include 'COMMON.INTERACT'
1250 include 'COMMON.LOCAL'
1251 include 'COMMON.NAMES'
1252 include 'COMMON.CHAIN'
1253 include 'COMMON.FFIELD'
1254 include 'COMMON.SBRIDGE'
1255 include 'COMMON.HEADER'
1256 include 'COMMON.CONTROL'
1257 include 'COMMON.DBASE'
1258 include 'COMMON.THREAD'
1259 include 'COMMON.TIME1'
1260 include 'COMMON.SETUP'
1261 C Read bridging residues.
1262 read (inp,*) ns,(iss(i),i=1,ns)
1264 if(me.eq.king.or..not.out1file)
1265 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1266 C Check whether the specified bridging residues are cystines.
1268 if (itype(iss(i)).ne.1) then
1269 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1270 & 'Do you REALLY think that the residue ',
1271 & restyp(itype(iss(i))),i,
1272 & ' can form a disulfide bridge?!!!'
1273 write (*,'(2a,i3,a)')
1274 & 'Do you REALLY think that the residue ',
1275 & restyp(itype(iss(i))),i,
1276 & ' can form a disulfide bridge?!!!'
1278 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1283 C Read preformed bridges.
1285 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1287 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1290 C Check if the residues involved in bridges are in the specified list of
1291 C bridging residues.
1294 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1295 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1296 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1297 & ' contains residues present in other pairs.'
1298 write (*,'(a,i3,a)') 'Disulfide pair',i,
1299 & ' contains residues present in other pairs.'
1301 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1307 if (ihpb(i).eq.iss(j)) goto 10
1309 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1312 if (jhpb(i).eq.iss(j)) goto 20
1314 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1320 ihpb(i)=ihpb(i)+nres
1321 jhpb(i)=jhpb(i)+nres
1327 c----------------------------------------------------------------------------
1328 subroutine read_x(kanal,*)
1329 implicit real*8 (a-h,o-z)
1330 include 'DIMENSIONS'
1331 include 'COMMON.GEO'
1332 include 'COMMON.VAR'
1333 include 'COMMON.CHAIN'
1334 include 'COMMON.IOUNITS'
1335 include 'COMMON.CONTROL'
1336 include 'COMMON.LOCAL'
1337 include 'COMMON.INTERACT'
1338 c Read coordinates from input
1340 read(kanal,'(8f10.5)',end=10,err=10)
1341 & ((c(l,k),l=1,3),k=1,nres),
1342 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1345 c(j,2*nres)=c(j,nres)
1347 call int_from_cart1(.false.)
1350 dc(j,i)=c(j,i+1)-c(j,i)
1351 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1355 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1357 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1358 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1366 c----------------------------------------------------------------------------
1367 subroutine read_threadbase
1368 implicit real*8 (a-h,o-z)
1369 include 'DIMENSIONS'
1370 include 'COMMON.IOUNITS'
1371 include 'COMMON.GEO'
1372 include 'COMMON.VAR'
1373 include 'COMMON.INTERACT'
1374 include 'COMMON.LOCAL'
1375 include 'COMMON.NAMES'
1376 include 'COMMON.CHAIN'
1377 include 'COMMON.FFIELD'
1378 include 'COMMON.SBRIDGE'
1379 include 'COMMON.HEADER'
1380 include 'COMMON.CONTROL'
1381 include 'COMMON.DBASE'
1382 include 'COMMON.THREAD'
1383 include 'COMMON.TIME1'
1384 C Read pattern database for threading.
1385 read (icbase,*) nseq
1387 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1388 & nres_base(2,i),nres_base(3,i)
1389 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1391 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1392 c & nres_base(2,i),nres_base(3,i)
1393 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1397 if (weidis.eq.0.0D0) weidis=0.1D0
1406 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1407 write (iout,'(a,i5)') 'nexcl: ',nexcl
1408 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1411 c------------------------------------------------------------------------------
1412 subroutine setup_var
1413 implicit real*8 (a-h,o-z)
1414 include 'DIMENSIONS'
1415 include 'COMMON.IOUNITS'
1416 include 'COMMON.GEO'
1417 include 'COMMON.VAR'
1418 include 'COMMON.INTERACT'
1419 include 'COMMON.LOCAL'
1420 include 'COMMON.NAMES'
1421 include 'COMMON.CHAIN'
1422 include 'COMMON.FFIELD'
1423 include 'COMMON.SBRIDGE'
1424 include 'COMMON.HEADER'
1425 include 'COMMON.CONTROL'
1426 include 'COMMON.DBASE'
1427 include 'COMMON.THREAD'
1428 include 'COMMON.TIME1'
1429 C Set up variable list.
1435 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1437 ialph(i,1)=nvar+nside
1441 if (indphi.gt.0) then
1443 else if (indback.gt.0) then
1448 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1451 c----------------------------------------------------------------------------
1452 subroutine gen_dist_constr
1453 C Generate CA distance constraints.
1454 implicit real*8 (a-h,o-z)
1455 include 'DIMENSIONS'
1456 include 'COMMON.IOUNITS'
1457 include 'COMMON.GEO'
1458 include 'COMMON.VAR'
1459 include 'COMMON.INTERACT'
1460 include 'COMMON.LOCAL'
1461 include 'COMMON.NAMES'
1462 include 'COMMON.CHAIN'
1463 include 'COMMON.FFIELD'
1464 include 'COMMON.SBRIDGE'
1465 include 'COMMON.HEADER'
1466 include 'COMMON.CONTROL'
1467 include 'COMMON.DBASE'
1468 include 'COMMON.THREAD'
1469 include 'COMMON.TIME1'
1470 dimension itype_pdb(maxres)
1471 common /pizda/ itype_pdb
1473 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1474 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1475 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1477 do i=nstart_sup,nstart_sup+nsup-1
1478 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1479 cd & ' seq_pdb', restyp(itype_pdb(i))
1480 do j=i+2,nstart_sup+nsup-1
1482 ihpb(nhpb)=i+nstart_seq-nstart_sup
1483 jhpb(nhpb)=j+nstart_seq-nstart_sup
1485 dhpb(nhpb)=dist(i,j)
1488 cd write (iout,'(a)') 'Distance constraints:'
1493 cd if (ii.gt.nres) then
1498 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1499 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1500 cd & dhpb(i),forcon(i)
1504 c----------------------------------------------------------------------------
1506 implicit real*8 (a-h,o-z)
1507 include 'DIMENSIONS'
1508 include 'COMMON.MAP'
1509 include 'COMMON.IOUNITS'
1510 character*3 angid(4) /'THE','PHI','ALP','OME'/
1511 character*80 mapcard,ucase
1513 read (inp,'(a)') mapcard
1514 mapcard=ucase(mapcard)
1515 if (index(mapcard,'PHI').gt.0) then
1517 else if (index(mapcard,'THE').gt.0) then
1519 else if (index(mapcard,'ALP').gt.0) then
1521 else if (index(mapcard,'OME').gt.0) then
1524 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1525 stop 'Error - illegal variable spec in MAP card.'
1527 call readi (mapcard,'RES1',res1(imap),0)
1528 call readi (mapcard,'RES2',res2(imap),0)
1529 if (res1(imap).eq.0) then
1530 res1(imap)=res2(imap)
1531 else if (res2(imap).eq.0) then
1532 res2(imap)=res1(imap)
1534 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1536 & 'Error - illegal definition of variable group in MAP.'
1537 stop 'Error - illegal definition of variable group in MAP.'
1539 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1540 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1541 call readi(mapcard,'NSTEP',nstep(imap),0)
1542 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1544 & 'Illegal boundary and/or step size specification in MAP.'
1545 stop 'Illegal boundary and/or step size specification in MAP.'
1550 c----------------------------------------------------------------------------
1552 implicit real*8 (a-h,o-z)
1553 include 'DIMENSIONS'
1554 include 'COMMON.IOUNITS'
1555 include 'COMMON.GEO'
1556 include 'COMMON.CSA'
1557 include 'COMMON.BANK'
1558 include 'COMMON.CONTROL'
1560 character*620 mcmcard
1561 call card_concat(mcmcard)
1563 call readi(mcmcard,'NCONF',nconf,50)
1564 call readi(mcmcard,'NADD',nadd,0)
1565 call readi(mcmcard,'JSTART',jstart,1)
1566 call readi(mcmcard,'JEND',jend,1)
1567 call readi(mcmcard,'NSTMAX',nstmax,500000)
1568 call readi(mcmcard,'N0',n0,1)
1569 call readi(mcmcard,'N1',n1,6)
1570 call readi(mcmcard,'N2',n2,4)
1571 call readi(mcmcard,'N3',n3,0)
1572 call readi(mcmcard,'N4',n4,0)
1573 call readi(mcmcard,'N5',n5,0)
1574 call readi(mcmcard,'N6',n6,10)
1575 call readi(mcmcard,'N7',n7,0)
1576 call readi(mcmcard,'N8',n8,0)
1577 call readi(mcmcard,'N9',n9,0)
1578 call readi(mcmcard,'N14',n14,0)
1579 call readi(mcmcard,'N15',n15,0)
1580 call readi(mcmcard,'N16',n16,0)
1581 call readi(mcmcard,'N17',n17,0)
1582 call readi(mcmcard,'N18',n18,0)
1584 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1586 call readi(mcmcard,'NDIFF',ndiff,2)
1587 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1588 call readi(mcmcard,'IS1',is1,1)
1589 call readi(mcmcard,'IS2',is2,8)
1590 call readi(mcmcard,'NRAN0',nran0,4)
1591 call readi(mcmcard,'NRAN1',nran1,2)
1592 call readi(mcmcard,'IRR',irr,1)
1593 call readi(mcmcard,'NSEED',nseed,20)
1594 call readi(mcmcard,'NTOTAL',ntotal,10000)
1595 call reada(mcmcard,'CUT1',cut1,2.0d0)
1596 call reada(mcmcard,'CUT2',cut2,5.0d0)
1597 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1598 call readi(mcmcard,'ICMAX',icmax,3)
1599 call readi(mcmcard,'IRESTART',irestart,0)
1600 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1603 call reada(mcmcard,'DELE',dele,20.0d0)
1604 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1605 call readi(mcmcard,'IREF',iref,0)
1606 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1607 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1608 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1609 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1610 write (iout,*) "NCONF_IN",nconf_in
1613 c----------------------------------------------------------------------------
1614 cfmc subroutine mcmfread
1615 cfmc implicit real*8 (a-h,o-z)
1616 cfmc include 'DIMENSIONS'
1617 cfmc include 'COMMON.MCMF'
1618 cfmc include 'COMMON.IOUNITS'
1619 cfmc include 'COMMON.GEO'
1620 cfmc character*80 ucase
1621 cfmc character*620 mcmcard
1622 cfmc call card_concat(mcmcard)
1624 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1625 cfmc write(iout,*)'MAXRANT=',maxrant
1626 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1627 cfmc write(iout,*)'MAXFAM=',maxfam
1628 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1629 cfmc write(iout,*)'NNET1=',nnet1
1630 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1631 cfmc write(iout,*)'NNET2=',nnet2
1632 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1633 cfmc write(iout,*)'NNET3=',nnet3
1634 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1635 cfmc write(iout,*)'ILASTT=',ilastt
1636 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1637 cfmc write(iout,*)'MAXSTR=',maxstr
1638 cfmc maxstr_f=maxstr/maxfam
1639 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1640 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1641 cfmc write(iout,*)'NMCMF=',nmcmf
1642 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1643 cfmc write(iout,*)'IFOCUS=',ifocus
1644 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1645 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1646 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1647 cfmc write(iout,*)'INTPRT=',intprt
1648 cfmc call readi(mcmcard,'IPRT',iprt,100)
1649 cfmc write(iout,*)'IPRT=',iprt
1650 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1651 cfmc write(iout,*)'IMAXTR=',imaxtr
1652 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1653 cfmc write(iout,*)'MAXEVEN=',maxeven
1654 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1655 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1656 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1657 cfmc write(iout,*)'INIMIN=',inimin
1658 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1659 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1660 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1661 cfmc write(iout,*)'NTHREAD=',nthread
1662 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1663 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1664 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1665 cfmc write(iout,*)'MAXPERT=',maxpert
1666 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1667 cfmc write(iout,*)'IRMSD=',irmsd
1668 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1669 cfmc write(iout,*)'DENEMIN=',denemin
1670 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1671 cfmc write(iout,*)'RCUT1S=',rcut1s
1672 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1673 cfmc write(iout,*)'RCUT1E=',rcut1e
1674 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1675 cfmc write(iout,*)'RCUT2S=',rcut2s
1676 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1677 cfmc write(iout,*)'RCUT2E=',rcut2e
1678 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1679 cfmc write(iout,*)'DPERT1=',d_pert1
1680 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1681 cfmc write(iout,*)'DPERT1A=',d_pert1a
1682 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1683 cfmc write(iout,*)'DPERT2=',d_pert2
1684 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1685 cfmc write(iout,*)'DPERT2A=',d_pert2a
1686 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1687 cfmc write(iout,*)'DPERT2B=',d_pert2b
1688 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1689 cfmc write(iout,*)'DPERT2C=',d_pert2c
1690 cfmc d_pert1=deg2rad*d_pert1
1691 cfmc d_pert1a=deg2rad*d_pert1a
1692 cfmc d_pert2=deg2rad*d_pert2
1693 cfmc d_pert2a=deg2rad*d_pert2a
1694 cfmc d_pert2b=deg2rad*d_pert2b
1695 cfmc d_pert2c=deg2rad*d_pert2c
1696 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1697 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1698 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1699 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1700 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1701 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1702 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1703 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1704 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1705 cfmc write(iout,*)'RCUTINI=',rcutini
1706 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1707 cfmc write(iout,*)'GRAT=',grat
1708 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1709 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1713 c----------------------------------------------------------------------------
1715 implicit real*8 (a-h,o-z)
1716 include 'DIMENSIONS'
1717 include 'COMMON.MCM'
1718 include 'COMMON.MCE'
1719 include 'COMMON.IOUNITS'
1721 character*320 mcmcard
1722 call card_concat(mcmcard)
1723 call readi(mcmcard,'MAXACC',maxacc,100)
1724 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1725 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1726 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1727 call readi(mcmcard,'MAXREPM',maxrepm,200)
1728 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1729 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1730 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1731 call reada(mcmcard,'E_UP',e_up,5.0D0)
1732 call reada(mcmcard,'DELTE',delte,0.1D0)
1733 call readi(mcmcard,'NSWEEP',nsweep,5)
1734 call readi(mcmcard,'NSTEPH',nsteph,0)
1735 call readi(mcmcard,'NSTEPC',nstepc,0)
1736 call reada(mcmcard,'TMIN',tmin,298.0D0)
1737 call reada(mcmcard,'TMAX',tmax,298.0D0)
1738 call readi(mcmcard,'NWINDOW',nwindow,0)
1739 call readi(mcmcard,'PRINT_MC',print_mc,0)
1740 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1741 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1742 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1743 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1744 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1745 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1746 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1747 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1748 if (nwindow.gt.0) then
1749 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1751 winlen(i)=winend(i)-winstart(i)+1
1754 if (tmax.lt.tmin) tmax=tmin
1755 if (tmax.eq.tmin) then
1759 if (nstepc.gt.0 .and. nsteph.gt.0) then
1760 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1761 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1763 C Probabilities of different move types
1764 sumpro_type(0)=0.0D0
1765 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1766 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1767 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1768 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1769 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1770 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1771 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1773 print *,'i',i,' sumprotype',sumpro_type(i)
1774 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1775 print *,'i',i,' sumprotype',sumpro_type(i)
1779 c----------------------------------------------------------------------------
1780 subroutine read_minim
1781 implicit real*8 (a-h,o-z)
1782 include 'DIMENSIONS'
1783 include 'COMMON.MINIM'
1784 include 'COMMON.IOUNITS'
1786 character*320 minimcard
1787 call card_concat(minimcard)
1788 call readi(minimcard,'MAXMIN',maxmin,2000)
1789 call readi(minimcard,'MAXFUN',maxfun,5000)
1790 call readi(minimcard,'MINMIN',minmin,maxmin)
1791 call readi(minimcard,'MINFUN',minfun,maxmin)
1792 call reada(minimcard,'TOLF',tolf,1.0D-2)
1793 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1794 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1795 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1796 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1797 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1798 & 'Options in energy minimization:'
1799 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1800 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1801 & 'MinMin:',MinMin,' MinFun:',MinFun,
1802 & ' TolF:',TolF,' RTolF:',RTolF
1805 c----------------------------------------------------------------------------
1806 subroutine read_angles(kanal,*)
1807 implicit real*8 (a-h,o-z)
1808 include 'DIMENSIONS'
1809 include 'COMMON.GEO'
1810 include 'COMMON.VAR'
1811 include 'COMMON.CHAIN'
1812 include 'COMMON.IOUNITS'
1813 include 'COMMON.CONTROL'
1814 c Read angles from input
1816 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1817 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1818 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1819 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1822 c 9/7/01 avoid 180 deg valence angle
1823 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1825 theta(i)=deg2rad*theta(i)
1826 phi(i)=deg2rad*phi(i)
1827 alph(i)=deg2rad*alph(i)
1828 omeg(i)=deg2rad*omeg(i)
1833 c----------------------------------------------------------------------------
1834 subroutine reada(rekord,lancuch,wartosc,default)
1836 character*(*) rekord,lancuch
1837 double precision wartosc,default
1840 iread=index(rekord,lancuch)
1841 if (iread.eq.0) then
1845 iread=iread+ilen(lancuch)+1
1846 read (rekord(iread:),*,err=10,end=10) wartosc
1851 c----------------------------------------------------------------------------
1852 subroutine readi(rekord,lancuch,wartosc,default)
1854 character*(*) rekord,lancuch
1855 integer wartosc,default
1858 iread=index(rekord,lancuch)
1859 if (iread.eq.0) then
1863 iread=iread+ilen(lancuch)+1
1864 read (rekord(iread:),*,err=10,end=10) wartosc
1869 c----------------------------------------------------------------------------
1870 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1873 integer tablica(dim),default
1874 character*(*) rekord,lancuch
1881 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1882 if (iread.eq.0) return
1883 iread=iread+ilen(lancuch)+1
1884 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1887 c----------------------------------------------------------------------------
1888 subroutine multreada(rekord,lancuch,tablica,dim,default)
1891 double precision tablica(dim),default
1892 character*(*) rekord,lancuch
1899 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1900 if (iread.eq.0) return
1901 iread=iread+ilen(lancuch)+1
1902 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1905 c----------------------------------------------------------------------------
1906 subroutine openunits
1907 implicit real*8 (a-h,o-z)
1908 include 'DIMENSIONS'
1911 character*16 form,nodename
1914 include 'COMMON.SETUP'
1915 include 'COMMON.IOUNITS'
1917 include 'COMMON.CONTROL'
1918 integer lenpre,lenpot,ilen,lentmp
1920 character*3 out1file_text,ucase
1923 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1924 call getenv_loc("PREFIX",prefix)
1926 call getenv_loc("POT",pot)
1927 call getenv_loc("DIRTMP",tmpdir)
1928 call getenv_loc("CURDIR",curdir)
1929 call getenv_loc("OUT1FILE",out1file_text)
1930 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1931 out1file_text=ucase(out1file_text)
1932 if (out1file_text(1:1).eq."Y") then
1935 out1file=fg_rank.gt.0
1940 if (lentmp.gt.0) then
1941 write (*,'(80(1h!))')
1942 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1943 write (*,'(80(1h!))')
1944 write (*,*)"All output files will be on node /tmp directory."
1946 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1947 if (me.eq.king) then
1948 write (*,*) "The master node is ",nodename
1949 else if (fg_rank.eq.0) then
1950 write (*,*) "I am the CG slave node ",nodename
1952 write (*,*) "I am the FG slave node ",nodename
1955 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1956 lenpre = lentmp+lenpre+1
1958 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1959 C Get the names and open the input files
1960 #if defined(WINIFL) || defined(WINPGI)
1961 open(1,file=pref_orig(:ilen(pref_orig))//
1962 & '.inp',status='old',readonly,shared)
1963 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1964 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1965 C Get parameter filenames and open the parameter files.
1966 call getenv_loc('BONDPAR',bondname)
1967 open (ibond,file=bondname,status='old',readonly,shared)
1968 call getenv_loc('THETPAR',thetname)
1969 open (ithep,file=thetname,status='old',readonly,shared)
1971 call getenv_loc('THETPARPDB',thetname_pdb)
1972 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1974 call getenv_loc('ROTPAR',rotname)
1975 open (irotam,file=rotname,status='old',readonly,shared)
1977 call getenv_loc('ROTPARPDB',rotname_pdb)
1978 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1980 call getenv_loc('TORPAR',torname)
1981 open (itorp,file=torname,status='old',readonly,shared)
1982 call getenv_loc('TORDPAR',tordname)
1983 open (itordp,file=tordname,status='old',readonly,shared)
1984 call getenv_loc('FOURIER',fouriername)
1985 open (ifourier,file=fouriername,status='old',readonly,shared)
1986 call getenv_loc('ELEPAR',elename)
1987 open (ielep,file=elename,status='old',readonly,shared)
1988 call getenv_loc('SIDEPAR',sidename)
1989 open (isidep,file=sidename,status='old',readonly,shared)
1990 #elif (defined CRAY) || (defined AIX)
1991 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1993 c print *,"Processor",myrank," opened file 1"
1994 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1995 c print *,"Processor",myrank," opened file 9"
1996 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1997 C Get parameter filenames and open the parameter files.
1998 call getenv_loc('BONDPAR',bondname)
1999 open (ibond,file=bondname,status='old',action='read')
2000 c print *,"Processor",myrank," opened file IBOND"
2001 call getenv_loc('THETPAR',thetname)
2002 open (ithep,file=thetname,status='old',action='read')
2003 c print *,"Processor",myrank," opened file ITHEP"
2005 call getenv_loc('THETPARPDB',thetname_pdb)
2006 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2008 call getenv_loc('ROTPAR',rotname)
2009 open (irotam,file=rotname,status='old',action='read')
2010 c print *,"Processor",myrank," opened file IROTAM"
2012 call getenv_loc('ROTPARPDB',rotname_pdb)
2013 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2015 call getenv_loc('TORPAR',torname)
2016 open (itorp,file=torname,status='old',action='read')
2017 c print *,"Processor",myrank," opened file ITORP"
2018 call getenv_loc('TORDPAR',tordname)
2019 open (itordp,file=tordname,status='old',action='read')
2020 c print *,"Processor",myrank," opened file ITORDP"
2021 call getenv_loc('SCCORPAR',sccorname)
2022 open (isccor,file=sccorname,status='old',action='read')
2023 c print *,"Processor",myrank," opened file ISCCOR"
2024 call getenv_loc('FOURIER',fouriername)
2025 open (ifourier,file=fouriername,status='old',action='read')
2026 c print *,"Processor",myrank," opened file IFOURIER"
2027 call getenv_loc('ELEPAR',elename)
2028 open (ielep,file=elename,status='old',action='read')
2029 c print *,"Processor",myrank," opened file IELEP"
2030 call getenv_loc('SIDEPAR',sidename)
2031 open (isidep,file=sidename,status='old',action='read')
2032 c print *,"Processor",myrank," opened file ISIDEP"
2033 c print *,"Processor",myrank," opened parameter files"
2035 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2036 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2037 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2038 C Get parameter filenames and open the parameter files.
2039 call getenv_loc('BONDPAR',bondname)
2040 open (ibond,file=bondname,status='old')
2041 call getenv_loc('THETPAR',thetname)
2042 open (ithep,file=thetname,status='old')
2044 call getenv_loc('THETPARPDB',thetname_pdb)
2045 open (ithep_pdb,file=thetname_pdb,status='old')
2047 call getenv_loc('ROTPAR',rotname)
2048 open (irotam,file=rotname,status='old')
2050 call getenv_loc('ROTPARPDB',rotname_pdb)
2051 open (irotam_pdb,file=rotname_pdb,status='old')
2053 call getenv_loc('TORPAR',torname)
2054 open (itorp,file=torname,status='old')
2055 call getenv_loc('TORDPAR',tordname)
2056 open (itordp,file=tordname,status='old')
2057 call getenv_loc('SCCORPAR',sccorname)
2058 open (isccor,file=sccorname,status='old')
2059 call getenv_loc('FOURIER',fouriername)
2060 open (ifourier,file=fouriername,status='old')
2061 call getenv_loc('ELEPAR',elename)
2062 open (ielep,file=elename,status='old')
2063 call getenv_loc('SIDEPAR',sidename)
2064 open (isidep,file=sidename,status='old')
2066 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2068 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2069 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2070 C Get parameter filenames and open the parameter files.
2071 call getenv_loc('BONDPAR',bondname)
2072 open (ibond,file=bondname,status='old',readonly)
2073 call getenv_loc('THETPAR',thetname)
2074 open (ithep,file=thetname,status='old',readonly)
2075 call getenv_loc('ROTPAR',rotname)
2076 open (irotam,file=rotname,status='old',readonly)
2077 call getenv_loc('TORPAR',torname)
2078 open (itorp,file=torname,status='old',readonly)
2079 call getenv_loc('TORDPAR',tordname)
2080 open (itordp,file=tordname,status='old',readonly)
2081 call getenv_loc('SCCORPAR',sccorname)
2082 open (isccor,file=sccorname,status='old',readonly)
2084 call getenv_loc('THETPARPDB',thetname_pdb)
2085 c print *,"thetname_pdb ",thetname_pdb
2086 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2087 c print *,ithep_pdb," opened"
2089 call getenv_loc('FOURIER',fouriername)
2090 open (ifourier,file=fouriername,status='old',readonly)
2091 call getenv_loc('ELEPAR',elename)
2092 open (ielep,file=elename,status='old',readonly)
2093 call getenv_loc('SIDEPAR',sidename)
2094 open (isidep,file=sidename,status='old',readonly)
2095 call getenv_loc('LIPTRANPAR',liptranname)
2096 open (iliptranpar,file=liptranname,status='old',action='read')
2098 call getenv_loc('ROTPARPDB',rotname_pdb)
2099 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2104 C 8/9/01 In the newest version SCp interaction constants are read from a file
2105 C Use -DOLDSCP to use hard-coded constants instead.
2107 call getenv_loc('SCPPAR',scpname)
2108 #if defined(WINIFL) || defined(WINPGI)
2109 open (iscpp,file=scpname,status='old',readonly,shared)
2110 #elif (defined CRAY) || (defined AIX)
2111 open (iscpp,file=scpname,status='old',action='read')
2113 open (iscpp,file=scpname,status='old')
2115 open (iscpp,file=scpname,status='old',readonly)
2118 call getenv_loc('PATTERN',patname)
2119 #if defined(WINIFL) || defined(WINPGI)
2120 open (icbase,file=patname,status='old',readonly,shared)
2121 #elif (defined CRAY) || (defined AIX)
2122 open (icbase,file=patname,status='old',action='read')
2124 open (icbase,file=patname,status='old')
2126 open (icbase,file=patname,status='old',readonly)
2129 C Open output file only for CG processes
2130 c print *,"Processor",myrank," fg_rank",fg_rank
2131 if (fg_rank.eq.0) then
2133 if (nodes.eq.1) then
2136 npos = dlog10(dfloat(nodes-1))+1
2138 if (npos.lt.3) npos=3
2139 write (liczba,'(i1)') npos
2140 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2142 write (liczba,form) me
2143 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2144 & liczba(:ilen(liczba))
2145 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2147 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2149 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2150 & liczba(:ilen(liczba))//'.mol2'
2151 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2152 & liczba(:ilen(liczba))//'.stat'
2154 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2155 & //liczba(:ilen(liczba))//'.stat')
2156 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2159 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2160 & liczba(:ilen(liczba))//'.const'
2165 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2166 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2167 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2168 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2169 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2171 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2173 rest2name=prefix(:ilen(prefix))//'.rst'
2175 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2178 #if defined(AIX) || defined(PGI)
2179 if (me.eq.king .or. .not. out1file)
2180 & open(iout,file=outname,status='unknown')
2182 if (fg_rank.gt.0) then
2183 write (liczba,'(i3.3)') myrank/nfgtasks
2184 write (ll,'(bz,i3.3)') fg_rank
2185 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2190 open(igeom,file=intname,status='unknown',position='append')
2191 open(ipdb,file=pdbname,status='unknown')
2192 open(imol2,file=mol2name,status='unknown')
2193 open(istat,file=statname,status='unknown',position='append')
2195 c1out open(iout,file=outname,status='unknown')
2198 if (me.eq.king .or. .not.out1file)
2199 & open(iout,file=outname,status='unknown')
2201 if (fg_rank.gt.0) then
2202 write (liczba,'(i3.3)') myrank/nfgtasks
2203 write (ll,'(bz,i3.3)') fg_rank
2204 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2209 open(igeom,file=intname,status='unknown',access='append')
2210 open(ipdb,file=pdbname,status='unknown')
2211 open(imol2,file=mol2name,status='unknown')
2212 open(istat,file=statname,status='unknown',access='append')
2214 c1out open(iout,file=outname,status='unknown')
2217 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2218 csa_seed=prefix(:lenpre)//'.CSA.seed'
2219 csa_history=prefix(:lenpre)//'.CSA.history'
2220 csa_bank=prefix(:lenpre)//'.CSA.bank'
2221 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2222 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2223 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2224 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2225 csa_int=prefix(:lenpre)//'.int'
2226 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2227 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2228 csa_in=prefix(:lenpre)//'.CSA.in'
2229 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2232 write (iout,'(80(1h-))')
2233 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2234 write (iout,'(80(1h-))')
2235 write (iout,*) "Input file : ",
2236 & pref_orig(:ilen(pref_orig))//'.inp'
2237 write (iout,*) "Output file : ",
2238 & outname(:ilen(outname))
2240 write (iout,*) "Sidechain potential file : ",
2241 & sidename(:ilen(sidename))
2243 write (iout,*) "SCp potential file : ",
2244 & scpname(:ilen(scpname))
2246 write (iout,*) "Electrostatic potential file : ",
2247 & elename(:ilen(elename))
2248 write (iout,*) "Cumulant coefficient file : ",
2249 & fouriername(:ilen(fouriername))
2250 write (iout,*) "Torsional parameter file : ",
2251 & torname(:ilen(torname))
2252 write (iout,*) "Double torsional parameter file : ",
2253 & tordname(:ilen(tordname))
2254 write (iout,*) "SCCOR parameter file : ",
2255 & sccorname(:ilen(sccorname))
2256 write (iout,*) "Bond & inertia constant file : ",
2257 & bondname(:ilen(bondname))
2258 write (iout,*) "Bending parameter file : ",
2259 & thetname(:ilen(thetname))
2260 write (iout,*) "Rotamer parameter file : ",
2261 & rotname(:ilen(rotname))
2262 write (iout,*) "Thetpdb parameter file : ",
2263 & thetname_pdb(:ilen(thetname_pdb))
2264 write (iout,*) "Threading database : ",
2265 & patname(:ilen(patname))
2267 &write (iout,*)" DIRTMP : ",
2269 write (iout,'(80(1h-))')
2273 c----------------------------------------------------------------------------
2274 subroutine card_concat(card)
2275 implicit real*8 (a-h,o-z)
2276 include 'DIMENSIONS'
2277 include 'COMMON.IOUNITS'
2279 character*80 karta,ucase
2281 read (inp,'(a)') karta
2284 do while (karta(80:80).eq.'&')
2285 card=card(:ilen(card)+1)//karta(:79)
2286 read (inp,'(a)') karta
2289 card=card(:ilen(card)+1)//karta
2292 c----------------------------------------------------------------------------------
2294 implicit real*8 (a-h,o-z)
2295 include 'DIMENSIONS'
2296 include 'COMMON.CHAIN'
2297 include 'COMMON.IOUNITS'
2299 include 'COMMON.CONTROL'
2300 open(irest2,file=rest2name,status='unknown')
2301 read(irest2,*) totT,EK,potE,totE,t_bath
2304 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2307 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2309 if(usampl.or.homol_nset.gt.1) then
2310 read (irest2,*) iset
2315 c---------------------------------------------------------------------------------
2316 subroutine read_fragments
2317 implicit real*8 (a-h,o-z)
2318 include 'DIMENSIONS'
2322 include 'COMMON.SETUP'
2323 include 'COMMON.CHAIN'
2324 include 'COMMON.IOUNITS'
2326 include 'COMMON.CONTROL'
2328 read(inp,*) nset,nfrag,npair,nfrag_back
2329 if(me.eq.king.or..not.out1file)
2330 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2331 & " nfrag_back",nfrag_back
2333 read(inp,*) mset(iset1)
2335 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2337 if(me.eq.king.or..not.out1file)
2338 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2339 & ifrag(2,i,iset1), qinfrag(i,iset1)
2342 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2344 if(me.eq.king.or..not.out1file)
2345 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2346 & ipair(2,i,iset1), qinpair(i,iset1)
2349 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2350 & wfrag_back(3,i,iset1),
2351 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2352 if(me.eq.king.or..not.out1file)
2353 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2354 & wfrag_back(2,i,iset1),
2355 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2356 & ifrag_back(2,i,iset1)
2361 C---------------------------------------------------------------------------
2362 subroutine read_afminp
2363 implicit real*8 (a-h,o-z)
2364 include 'DIMENSIONS'
2368 include 'COMMON.SETUP'
2369 include 'COMMON.CONTROL'
2370 include 'COMMON.CHAIN'
2371 include 'COMMON.IOUNITS'
2372 include 'COMMON.SBRIDGE'
2373 character*320 afmcard
2374 c print *, "wchodze"
2375 call card_concat(afmcard)
2376 call readi(afmcard,"BEG",afmbeg,0)
2377 call readi(afmcard,"END",afmend,0)
2378 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2379 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2380 print *,'FORCE=' ,forceAFMconst
2381 CCCC NOW PROPERTIES FOR AFM
2384 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2386 distafminit=dsqrt(distafminit)
2387 print *,'initdist',distafminit
2391 c-------------------------------------------------------------------------------
2392 subroutine read_dist_constr
2393 implicit real*8 (a-h,o-z)
2394 include 'DIMENSIONS'
2398 include 'COMMON.SETUP'
2399 include 'COMMON.CONTROL'
2400 include 'COMMON.CHAIN'
2401 include 'COMMON.IOUNITS'
2402 include 'COMMON.SBRIDGE'
2403 integer ifrag_(2,100),ipair_(2,100)
2404 double precision wfrag_(100),wpair_(100)
2405 character*500 controlcard
2406 c write (iout,*) "Calling read_dist_constr"
2407 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2409 call card_concat(controlcard)
2410 call readi(controlcard,"NFRAG",nfrag_,0)
2411 call readi(controlcard,"NPAIR",npair_,0)
2412 call readi(controlcard,"NDIST",ndist_,0)
2413 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2414 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2415 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2416 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2417 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2418 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2419 c write (iout,*) "IFRAG"
2421 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2423 c write (iout,*) "IPAIR"
2425 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2427 if (.not.refstr .and. nfrag.gt.0) then
2429 & "ERROR: no reference structure to compute distance restraints"
2431 & "Restraints must be specified explicitly (NDIST=number)"
2434 if (nfrag.lt.2 .and. npair.gt.0) then
2435 write (iout,*) "ERROR: Less than 2 fragments specified",
2436 & " but distance restraints between pairs requested"
2441 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2442 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2443 & ifrag_(2,i)=nstart_sup+nsup-1
2444 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2446 if (wfrag_(i).gt.0.0d0) then
2447 do j=ifrag_(1,i),ifrag_(2,i)-1
2448 do k=j+1,ifrag_(2,i)
2449 c write (iout,*) "j",j," k",k
2451 if (constr_dist.eq.1) then
2456 forcon(nhpb)=wfrag_(i)
2457 else if (constr_dist.eq.2) then
2458 if (ddjk.le.dist_cut) then
2463 forcon(nhpb)=wfrag_(i)
2470 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2473 if (.not.out1file .or. me.eq.king)
2474 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2475 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2477 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2478 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2485 if (wpair_(i).gt.0.0d0) then
2493 do j=ifrag_(1,ii),ifrag_(2,ii)
2494 do k=ifrag_(1,jj),ifrag_(2,jj)
2498 forcon(nhpb)=wpair_(i)
2499 dhpb(nhpb)=dist(j,k)
2501 if (.not.out1file .or. me.eq.king)
2502 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2503 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2505 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2506 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2513 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2514 if (forcon(nhpb+1).gt.0.0d0) then
2516 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2518 if (.not.out1file .or. me.eq.king)
2519 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2520 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2522 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2523 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2530 c-------------------------------------------------------------------------------
2532 subroutine read_constr_homology
2534 include 'DIMENSIONS'
2538 include 'COMMON.SETUP'
2539 include 'COMMON.CONTROL'
2540 include 'COMMON.CHAIN'
2541 include 'COMMON.IOUNITS'
2543 include 'COMMON.GEO'
2544 include 'COMMON.INTERACT'
2546 c For new homol impl
2548 include 'COMMON.VAR'
2551 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2553 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2554 c & sigma_odl_temp(maxres,maxres,max_template)
2556 character*24 model_ki_dist, model_ki_angle
2557 character*500 controlcard
2558 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2559 logical lprn /.true./
2563 c FP - Nov. 2014 Temporary specifications for new vars
2565 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
2566 double precision, dimension (max_template,maxres) :: rescore
2567 double precision, dimension (max_template,maxres) :: rescore2
2568 character*24 pdbfile,tpl_k_rescore
2569 c -----------------------------------------------------------------
2570 c Reading multiple PDB ref structures and calculation of retraints
2571 c not using pre-computed ones stored in files model_ki_{dist,angle}
2573 c -----------------------------------------------------------------
2576 c Alternative: reading from input
2577 call card_concat(controlcard)
2578 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2579 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2580 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2581 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2582 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2583 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2584 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2585 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2586 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2587 if(.not.read2sigma.and.start_from_model) then
2588 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2589 start_from_model=.false.
2591 if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
2592 if(start_from_model .and. rest) then
2593 write(iout,*) 'START_FROM_MODELS is OFF'
2594 write(iout,*) 'remove restart keyword from input'
2596 if (homol_nset.gt.1)then
2597 call card_concat(controlcard)
2598 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2599 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2600 write(iout,*) "iset homology_weight "
2602 write(iout,*) i,waga_homology(i)
2605 iset=mod(kolor,homol_nset)+1
2608 waga_homology(1)=1.0
2611 cd write (iout,*) "nnt",nnt," nct",nct
2623 write(iout,*) 'nnt=',nnt,'nct=',nct
2626 do k=1,constr_homology
2639 do k=1,constr_homology
2641 read(inp,'(a)') pdbfile
2642 c Next stament causes error upon compilation (?)
2643 c if(me.eq.king.or. .not. out1file)
2644 c write (iout,'(2a)') 'PDB data will be read from file ',
2645 c & pdbfile(:ilen(pdbfile))
2646 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2647 & pdbfile(:ilen(pdbfile))
2648 open(ipdbin,file=pdbfile,status='old',err=33)
2650 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2651 & pdbfile(:ilen(pdbfile))
2654 c print *,'Begin reading pdb data'
2656 c Files containing res sim or local scores (former containing sigmas)
2659 write(kic2,'(bz,i2.2)') k
2661 tpl_k_rescore="template"//kic2//".sco"
2664 if (read2sigma) then
2665 call readpdb_template(k)
2670 c Distance restraints
2673 C Copy the coordinates from reference coordinates (?)
2677 c write (iout,*) "c(",j,i,") =",c(j,i)
2681 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2683 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2684 open (ientin,file=tpl_k_rescore,status='old')
2685 if (nnt.gt.1) rescore(k,1)=0.0d0
2686 do irec=nnt,maxdim ! loop for reading res sim
2687 if (read2sigma) then
2688 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2691 idomain(k,i_tmp)=idomain_tmp
2692 rescore(k,i_tmp)=rescore_tmp
2693 rescore2(k,i_tmp)=rescore2_tmp
2696 read (ientin,*,end=1401) rescore_tmp
2698 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2699 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2700 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2705 if (waga_dist.ne.0.0d0) then
2713 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2714 c write (iout,*) k,i,j,distal,dist2_cut
2716 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2717 & .and. distal.le.dist2_cut ) then
2723 c write (iout,*) "k",k
2724 c write (iout,*) "i",i," j",j," constr_homology",
2729 if (read2sigma) then
2732 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2734 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2735 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2736 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2738 if (odl(k,ii).le.dist_cut) then
2739 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2742 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2743 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2745 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2746 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2750 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2753 l_homo(k,ii)=.false.
2760 c Theta, dihedral and SC retraints
2762 if (waga_angle.gt.0.0d0) then
2763 c open (ientin,file=tpl_k_sigma_dih,status='old')
2764 c do irec=1,maxres-3 ! loop for reading sigma_dih
2765 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2766 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2767 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2768 c & sigma_dih(k,i+nnt-1)
2773 if (idomain(k,i).eq.0) then
2777 dih(k,i)=phiref(i) ! right?
2778 c read (ientin,*) sigma_dih(k,i) ! original variant
2779 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2780 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2781 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2782 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2783 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2785 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2786 & rescore(k,i-2)+rescore(k,i-3))/4.0
2787 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2788 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2789 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2790 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2791 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2792 c Instead of res sim other local measure of b/b str reliability possible
2793 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2794 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2799 if (waga_theta.gt.0.0d0) then
2800 c open (ientin,file=tpl_k_sigma_theta,status='old')
2801 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2802 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2803 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2804 c & sigma_theta(k,i+nnt-1)
2809 do i = nnt+2,nct ! right? without parallel.
2810 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2811 c do i=ithet_start,ithet_end ! with FG parallel.
2812 if (idomain(k,i).eq.0) then
2813 sigma_theta(k,i)=0.0
2816 thetatpl(k,i)=thetaref(i)
2817 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2818 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2819 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2820 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2821 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2822 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2823 & rescore(k,i-2))/3.0
2824 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2825 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2827 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2828 c rescore(k,i-2) ! right expression ?
2829 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2834 if (waga_d.gt.0.0d0) then
2835 c open (ientin,file=tpl_k_sigma_d,status='old')
2836 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2837 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2838 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2839 c & sigma_d(k,i+nnt-1)
2843 do i = nnt,nct ! right? without parallel.
2844 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2845 c do i=loc_start,loc_end ! with FG parallel.
2846 if (itype(i).eq.10) cycle
2847 if (idomain(k,i).eq.0 ) then
2854 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2855 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2856 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2857 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
2858 sigma_d(k,i)=rescore(k,i) ! right expression ?
2859 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2861 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
2862 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2863 c read (ientin,*) sigma_d(k,i) ! 1st variant
2864 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
2870 c remove distance restraints not used in any model from the list
2871 c shift data in all arrays
2873 if (waga_dist.ne.0.0d0) then
2878 if (ii_in_use(ii).eq.0) then
2880 ires_homo(ki)=ires_homo(ki+1)
2881 jres_homo(ki)=jres_homo(ki+1)
2882 ii_in_use(ki)=ii_in_use(ki+1)
2883 do k=1,constr_homology
2884 odl(k,ki)=odl(k,ki+1)
2885 sigma_odl(k,ki)=sigma_odl(k,ki+1)
2886 l_homo(k,ki)=l_homo(k,ki+1)
2895 if (constr_homology.gt.0) call homology_partition
2896 if (constr_homology.gt.0) call init_int_table
2897 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
2898 cd & "lim_xx=",lim_xx
2899 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2900 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2904 if (.not.lprn) return
2905 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2906 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2907 write (iout,*) "Distance restraints from templates"
2909 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
2910 & ii,ires_homo(ii),jres_homo(ii),
2911 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
2912 & ki=1,constr_homology)
2914 write (iout,*) "Dihedral angle restraints from templates"
2916 write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
2917 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2919 write (iout,*) "Virtual-bond angle restraints from templates"
2920 do i=nnt+2,lim_theta
2921 write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
2922 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2924 write (iout,*) "SC restraints from templates"
2926 write(iout,'(i5,100(4f8.2,4x))') i,
2927 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
2928 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2931 c -----------------------------------------------------------------
2934 c----------------------------------------------------------------------
2937 subroutine flush(iu)
2942 subroutine flush(iu)
2947 c------------------------------------------------------------------------------
2948 subroutine copy_to_tmp(source)
2949 include "DIMENSIONS"
2950 include "COMMON.IOUNITS"
2951 character*(*) source
2952 character* 256 tmpfile
2956 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2957 inquire(file=tmpfile,exist=ex)
2959 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2960 & " to temporary directory..."
2961 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2962 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2966 c------------------------------------------------------------------------------
2967 subroutine move_from_tmp(source)
2968 include "DIMENSIONS"
2969 include "COMMON.IOUNITS"
2970 character*(*) source
2973 write (*,*) "Moving ",source(:ilen(source)),
2974 & " from temporary directory to working directory"
2975 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2976 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2979 c------------------------------------------------------------------------------
2980 subroutine random_init(seed)
2982 C Initialize random number generator
2984 implicit real*8 (a-h,o-z)
2985 include 'DIMENSIONS'
2991 logical OKRandom, prng_restart
2993 integer iseed_array(4)
2995 include 'COMMON.IOUNITS'
2996 include 'COMMON.TIME1'
2997 include 'COMMON.THREAD'
2998 include 'COMMON.SBRIDGE'
2999 include 'COMMON.CONTROL'
3000 include 'COMMON.MCM'
3001 include 'COMMON.MAP'
3002 include 'COMMON.HEADER'
3003 include 'COMMON.CSA'
3004 include 'COMMON.CHAIN'
3005 include 'COMMON.MUCA'
3007 include 'COMMON.FFIELD'
3008 include 'COMMON.SETUP'
3009 iseed=-dint(dabs(seed))
3010 if (iseed.eq.0) then
3011 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3012 & 'Random seed undefined. The program will stop.'
3013 write (*,'(/80(1h*)/20x,a/80(1h*))')
3014 & 'Random seed undefined. The program will stop.'
3016 call mpi_finalize(mpi_comm_world,ierr)
3018 stop 'Bad random seed.'
3021 if (fg_rank.eq.0) then
3025 if(me.eq.king .or. .not. out1file)
3026 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3027 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3028 OKRandom = prng_restart(me,iseedi8)
3031 tmp=65536.0d0**(4-i)
3032 iseed_array(i) = dint(seed/tmp)
3033 seed=seed-iseed_array(i)*tmp
3036 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3037 & (iseed_array(i),i=1,4)
3038 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3039 & (iseed_array(i),i=1,4)
3040 OKRandom = prng_restart(me,iseed_array)
3043 c r1 = prng_next(me)
3044 r1=ran_number(0.0D0,1.0D0)
3045 if(me.eq.king .or. .not. out1file)
3046 & write (iout,*) 'ran_num',r1
3047 if (r1.lt.0.0d0) OKRandom=.false.
3049 if (.not.OKRandom) then
3050 write (iout,*) 'PRNG IS NOT WORKING!!!'
3051 print *,'PRNG IS NOT WORKING!!!'
3054 call mpi_abort(mpi_comm_world,error_msg,ierr)
3057 write (iout,*) 'too many processors for parallel prng'
3058 write (*,*) 'too many processors for parallel prng'
3066 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)