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 c write(iout,*) "bordliptop=",bordliptop
245 c write(iout,*) "bordlipbot=",bordlipbot
246 c write(iout,*) "bufliptop=",bufliptop
247 c 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 c 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'
1785 include 'COMMON.CONTROL'
1786 include 'COMMON.SETUP'
1788 character*320 minimcard
1789 call card_concat(minimcard)
1790 call readi(minimcard,'MAXMIN',maxmin,2000)
1791 call readi(minimcard,'MAXFUN',maxfun,5000)
1792 call readi(minimcard,'MINMIN',minmin,maxmin)
1793 call readi(minimcard,'MINFUN',minfun,maxmin)
1794 call reada(minimcard,'TOLF',tolf,1.0D-2)
1795 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1796 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1797 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1798 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1800 if (.not. out1file .or. me.eq.king) then
1802 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1803 & 'Options in energy minimization:'
1804 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1805 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1806 & 'MinMin:',MinMin,' MinFun:',MinFun,
1807 & ' TolF:',TolF,' RTolF:',RTolF
1813 c----------------------------------------------------------------------------
1814 subroutine read_angles(kanal,*)
1815 implicit real*8 (a-h,o-z)
1816 include 'DIMENSIONS'
1817 include 'COMMON.GEO'
1818 include 'COMMON.VAR'
1819 include 'COMMON.CHAIN'
1820 include 'COMMON.IOUNITS'
1821 include 'COMMON.CONTROL'
1822 c Read angles from input
1824 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1825 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1826 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1827 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1830 c 9/7/01 avoid 180 deg valence angle
1831 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1833 theta(i)=deg2rad*theta(i)
1834 phi(i)=deg2rad*phi(i)
1835 alph(i)=deg2rad*alph(i)
1836 omeg(i)=deg2rad*omeg(i)
1841 c----------------------------------------------------------------------------
1842 subroutine reada(rekord,lancuch,wartosc,default)
1844 character*(*) rekord,lancuch
1845 double precision wartosc,default
1848 iread=index(rekord,lancuch)
1849 if (iread.eq.0) then
1853 iread=iread+ilen(lancuch)+1
1854 read (rekord(iread:),*,err=10,end=10) wartosc
1859 c----------------------------------------------------------------------------
1860 subroutine readi(rekord,lancuch,wartosc,default)
1862 character*(*) rekord,lancuch
1863 integer wartosc,default
1866 iread=index(rekord,lancuch)
1867 if (iread.eq.0) then
1871 iread=iread+ilen(lancuch)+1
1872 read (rekord(iread:),*,err=10,end=10) wartosc
1877 c----------------------------------------------------------------------------
1878 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1881 integer tablica(dim),default
1882 character*(*) rekord,lancuch
1889 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1890 if (iread.eq.0) return
1891 iread=iread+ilen(lancuch)+1
1892 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1895 c----------------------------------------------------------------------------
1896 subroutine multreada(rekord,lancuch,tablica,dim,default)
1899 double precision tablica(dim),default
1900 character*(*) rekord,lancuch
1907 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1908 if (iread.eq.0) return
1909 iread=iread+ilen(lancuch)+1
1910 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1913 c----------------------------------------------------------------------------
1914 subroutine openunits
1915 implicit real*8 (a-h,o-z)
1916 include 'DIMENSIONS'
1919 character*16 form,nodename
1922 include 'COMMON.SETUP'
1923 include 'COMMON.IOUNITS'
1925 include 'COMMON.CONTROL'
1926 integer lenpre,lenpot,ilen,lentmp
1928 character*3 out1file_text,ucase
1931 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1932 call getenv_loc("PREFIX",prefix)
1934 call getenv_loc("POT",pot)
1935 call getenv_loc("DIRTMP",tmpdir)
1936 call getenv_loc("CURDIR",curdir)
1937 call getenv_loc("OUT1FILE",out1file_text)
1938 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1939 out1file_text=ucase(out1file_text)
1940 if (out1file_text(1:1).eq."Y") then
1943 out1file=fg_rank.gt.0
1948 if (lentmp.gt.0) then
1949 write (*,'(80(1h!))')
1950 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1951 write (*,'(80(1h!))')
1952 write (*,*)"All output files will be on node /tmp directory."
1954 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1955 if (me.eq.king) then
1956 write (*,*) "The master node is ",nodename
1957 else if (fg_rank.eq.0) then
1958 write (*,*) "I am the CG slave node ",nodename
1960 write (*,*) "I am the FG slave node ",nodename
1963 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1964 lenpre = lentmp+lenpre+1
1966 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1967 C Get the names and open the input files
1968 #if defined(WINIFL) || defined(WINPGI)
1969 open(1,file=pref_orig(:ilen(pref_orig))//
1970 & '.inp',status='old',readonly,shared)
1971 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1972 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1973 C Get parameter filenames and open the parameter files.
1974 call getenv_loc('BONDPAR',bondname)
1975 open (ibond,file=bondname,status='old',readonly,shared)
1976 call getenv_loc('THETPAR',thetname)
1977 open (ithep,file=thetname,status='old',readonly,shared)
1979 call getenv_loc('THETPARPDB',thetname_pdb)
1980 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1982 call getenv_loc('ROTPAR',rotname)
1983 open (irotam,file=rotname,status='old',readonly,shared)
1985 call getenv_loc('ROTPARPDB',rotname_pdb)
1986 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1988 call getenv_loc('TORPAR',torname)
1989 open (itorp,file=torname,status='old',readonly,shared)
1990 call getenv_loc('TORDPAR',tordname)
1991 open (itordp,file=tordname,status='old',readonly,shared)
1992 call getenv_loc('FOURIER',fouriername)
1993 open (ifourier,file=fouriername,status='old',readonly,shared)
1994 call getenv_loc('ELEPAR',elename)
1995 open (ielep,file=elename,status='old',readonly,shared)
1996 call getenv_loc('SIDEPAR',sidename)
1997 open (isidep,file=sidename,status='old',readonly,shared)
1998 #elif (defined CRAY) || (defined AIX)
1999 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2001 c print *,"Processor",myrank," opened file 1"
2002 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2003 c print *,"Processor",myrank," opened file 9"
2004 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2005 C Get parameter filenames and open the parameter files.
2006 call getenv_loc('BONDPAR',bondname)
2007 open (ibond,file=bondname,status='old',action='read')
2008 c print *,"Processor",myrank," opened file IBOND"
2009 call getenv_loc('THETPAR',thetname)
2010 open (ithep,file=thetname,status='old',action='read')
2011 c print *,"Processor",myrank," opened file ITHEP"
2013 call getenv_loc('THETPARPDB',thetname_pdb)
2014 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2016 call getenv_loc('ROTPAR',rotname)
2017 open (irotam,file=rotname,status='old',action='read')
2018 c print *,"Processor",myrank," opened file IROTAM"
2020 call getenv_loc('ROTPARPDB',rotname_pdb)
2021 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2023 call getenv_loc('TORPAR',torname)
2024 open (itorp,file=torname,status='old',action='read')
2025 c print *,"Processor",myrank," opened file ITORP"
2026 call getenv_loc('TORDPAR',tordname)
2027 open (itordp,file=tordname,status='old',action='read')
2028 c print *,"Processor",myrank," opened file ITORDP"
2029 call getenv_loc('SCCORPAR',sccorname)
2030 open (isccor,file=sccorname,status='old',action='read')
2031 c print *,"Processor",myrank," opened file ISCCOR"
2032 call getenv_loc('FOURIER',fouriername)
2033 open (ifourier,file=fouriername,status='old',action='read')
2034 c print *,"Processor",myrank," opened file IFOURIER"
2035 call getenv_loc('ELEPAR',elename)
2036 open (ielep,file=elename,status='old',action='read')
2037 c print *,"Processor",myrank," opened file IELEP"
2038 call getenv_loc('SIDEPAR',sidename)
2039 open (isidep,file=sidename,status='old',action='read')
2040 c print *,"Processor",myrank," opened file ISIDEP"
2041 c print *,"Processor",myrank," opened parameter files"
2043 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2044 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2045 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2046 C Get parameter filenames and open the parameter files.
2047 call getenv_loc('BONDPAR',bondname)
2048 open (ibond,file=bondname,status='old')
2049 call getenv_loc('THETPAR',thetname)
2050 open (ithep,file=thetname,status='old')
2052 call getenv_loc('THETPARPDB',thetname_pdb)
2053 open (ithep_pdb,file=thetname_pdb,status='old')
2055 call getenv_loc('ROTPAR',rotname)
2056 open (irotam,file=rotname,status='old')
2058 call getenv_loc('ROTPARPDB',rotname_pdb)
2059 open (irotam_pdb,file=rotname_pdb,status='old')
2061 call getenv_loc('TORPAR',torname)
2062 open (itorp,file=torname,status='old')
2063 call getenv_loc('TORDPAR',tordname)
2064 open (itordp,file=tordname,status='old')
2065 call getenv_loc('SCCORPAR',sccorname)
2066 open (isccor,file=sccorname,status='old')
2067 call getenv_loc('FOURIER',fouriername)
2068 open (ifourier,file=fouriername,status='old')
2069 call getenv_loc('ELEPAR',elename)
2070 open (ielep,file=elename,status='old')
2071 call getenv_loc('SIDEPAR',sidename)
2072 open (isidep,file=sidename,status='old')
2074 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2076 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2077 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2078 C Get parameter filenames and open the parameter files.
2079 call getenv_loc('BONDPAR',bondname)
2080 open (ibond,file=bondname,status='old',readonly)
2081 call getenv_loc('THETPAR',thetname)
2082 open (ithep,file=thetname,status='old',readonly)
2083 call getenv_loc('ROTPAR',rotname)
2084 open (irotam,file=rotname,status='old',readonly)
2085 call getenv_loc('TORPAR',torname)
2086 open (itorp,file=torname,status='old',readonly)
2087 call getenv_loc('TORDPAR',tordname)
2088 open (itordp,file=tordname,status='old',readonly)
2089 call getenv_loc('SCCORPAR',sccorname)
2090 open (isccor,file=sccorname,status='old',readonly)
2092 call getenv_loc('THETPARPDB',thetname_pdb)
2093 c print *,"thetname_pdb ",thetname_pdb
2094 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2095 c print *,ithep_pdb," opened"
2097 call getenv_loc('FOURIER',fouriername)
2098 open (ifourier,file=fouriername,status='old',readonly)
2099 call getenv_loc('ELEPAR',elename)
2100 open (ielep,file=elename,status='old',readonly)
2101 call getenv_loc('SIDEPAR',sidename)
2102 open (isidep,file=sidename,status='old',readonly)
2103 call getenv_loc('LIPTRANPAR',liptranname)
2104 open (iliptranpar,file=liptranname,status='old',action='read')
2106 call getenv_loc('ROTPARPDB',rotname_pdb)
2107 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2112 C 8/9/01 In the newest version SCp interaction constants are read from a file
2113 C Use -DOLDSCP to use hard-coded constants instead.
2115 call getenv_loc('SCPPAR',scpname)
2116 #if defined(WINIFL) || defined(WINPGI)
2117 open (iscpp,file=scpname,status='old',readonly,shared)
2118 #elif (defined CRAY) || (defined AIX)
2119 open (iscpp,file=scpname,status='old',action='read')
2121 open (iscpp,file=scpname,status='old')
2123 open (iscpp,file=scpname,status='old',readonly)
2126 call getenv_loc('PATTERN',patname)
2127 #if defined(WINIFL) || defined(WINPGI)
2128 open (icbase,file=patname,status='old',readonly,shared)
2129 #elif (defined CRAY) || (defined AIX)
2130 open (icbase,file=patname,status='old',action='read')
2132 open (icbase,file=patname,status='old')
2134 open (icbase,file=patname,status='old',readonly)
2137 C Open output file only for CG processes
2138 c print *,"Processor",myrank," fg_rank",fg_rank
2139 if (fg_rank.eq.0) then
2141 if (nodes.eq.1) then
2144 npos = dlog10(dfloat(nodes-1))+1
2146 if (npos.lt.3) npos=3
2147 write (liczba,'(i1)') npos
2148 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2150 write (liczba,form) me
2151 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2152 & liczba(:ilen(liczba))
2153 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2155 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2157 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2158 & liczba(:ilen(liczba))//'.mol2'
2159 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2160 & liczba(:ilen(liczba))//'.stat'
2162 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2163 & //liczba(:ilen(liczba))//'.stat')
2164 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2167 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2168 & liczba(:ilen(liczba))//'.const'
2173 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2174 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2175 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2176 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2177 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2179 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2181 rest2name=prefix(:ilen(prefix))//'.rst'
2183 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2186 #if defined(AIX) || defined(PGI)
2187 if (me.eq.king .or. .not. out1file)
2188 & open(iout,file=outname,status='unknown')
2190 if (fg_rank.gt.0) then
2191 write (liczba,'(i3.3)') myrank/nfgtasks
2192 write (ll,'(bz,i3.3)') fg_rank
2193 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2198 open(igeom,file=intname,status='unknown',position='append')
2199 open(ipdb,file=pdbname,status='unknown')
2200 open(imol2,file=mol2name,status='unknown')
2201 open(istat,file=statname,status='unknown',position='append')
2203 c1out open(iout,file=outname,status='unknown')
2206 if (me.eq.king .or. .not.out1file)
2207 & open(iout,file=outname,status='unknown')
2209 if (fg_rank.gt.0) then
2210 write (liczba,'(i3.3)') myrank/nfgtasks
2211 write (ll,'(bz,i3.3)') fg_rank
2212 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2217 open(igeom,file=intname,status='unknown',access='append')
2218 open(ipdb,file=pdbname,status='unknown')
2219 open(imol2,file=mol2name,status='unknown')
2220 open(istat,file=statname,status='unknown',access='append')
2222 c1out open(iout,file=outname,status='unknown')
2225 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2226 csa_seed=prefix(:lenpre)//'.CSA.seed'
2227 csa_history=prefix(:lenpre)//'.CSA.history'
2228 csa_bank=prefix(:lenpre)//'.CSA.bank'
2229 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2230 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2231 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2232 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2233 csa_int=prefix(:lenpre)//'.int'
2234 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2235 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2236 csa_in=prefix(:lenpre)//'.CSA.in'
2237 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2240 write (iout,'(80(1h-))')
2241 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2242 write (iout,'(80(1h-))')
2243 write (iout,*) "Input file : ",
2244 & pref_orig(:ilen(pref_orig))//'.inp'
2245 write (iout,*) "Output file : ",
2246 & outname(:ilen(outname))
2248 write (iout,*) "Sidechain potential file : ",
2249 & sidename(:ilen(sidename))
2251 write (iout,*) "SCp potential file : ",
2252 & scpname(:ilen(scpname))
2254 write (iout,*) "Electrostatic potential file : ",
2255 & elename(:ilen(elename))
2256 write (iout,*) "Cumulant coefficient file : ",
2257 & fouriername(:ilen(fouriername))
2258 write (iout,*) "Torsional parameter file : ",
2259 & torname(:ilen(torname))
2260 write (iout,*) "Double torsional parameter file : ",
2261 & tordname(:ilen(tordname))
2262 write (iout,*) "SCCOR parameter file : ",
2263 & sccorname(:ilen(sccorname))
2264 write (iout,*) "Bond & inertia constant file : ",
2265 & bondname(:ilen(bondname))
2266 write (iout,*) "Bending parameter file : ",
2267 & thetname(:ilen(thetname))
2268 write (iout,*) "Rotamer parameter file : ",
2269 & rotname(:ilen(rotname))
2270 write (iout,*) "Thetpdb parameter file : ",
2271 & thetname_pdb(:ilen(thetname_pdb))
2272 write (iout,*) "Threading database : ",
2273 & patname(:ilen(patname))
2275 &write (iout,*)" DIRTMP : ",
2277 write (iout,'(80(1h-))')
2281 c----------------------------------------------------------------------------
2282 subroutine card_concat(card)
2283 implicit real*8 (a-h,o-z)
2284 include 'DIMENSIONS'
2285 include 'COMMON.IOUNITS'
2287 character*80 karta,ucase
2289 read (inp,'(a)') karta
2292 do while (karta(80:80).eq.'&')
2293 card=card(:ilen(card)+1)//karta(:79)
2294 read (inp,'(a)') karta
2297 card=card(:ilen(card)+1)//karta
2300 c----------------------------------------------------------------------------------
2302 implicit real*8 (a-h,o-z)
2303 include 'DIMENSIONS'
2304 include 'COMMON.CHAIN'
2305 include 'COMMON.IOUNITS'
2307 include 'COMMON.CONTROL'
2308 open(irest2,file=rest2name,status='unknown')
2309 read(irest2,*) totT,EK,potE,totE,t_bath
2312 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2315 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2317 if(usampl.or.homol_nset.gt.1) then
2318 read (irest2,*) iset
2323 c---------------------------------------------------------------------------------
2324 subroutine read_fragments
2325 implicit real*8 (a-h,o-z)
2326 include 'DIMENSIONS'
2330 include 'COMMON.SETUP'
2331 include 'COMMON.CHAIN'
2332 include 'COMMON.IOUNITS'
2334 include 'COMMON.CONTROL'
2336 read(inp,*) nset,nfrag,npair,nfrag_back
2337 if(me.eq.king.or..not.out1file)
2338 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2339 & " nfrag_back",nfrag_back
2341 read(inp,*) mset(iset1)
2343 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2345 if(me.eq.king.or..not.out1file)
2346 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2347 & ifrag(2,i,iset1), qinfrag(i,iset1)
2350 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2352 if(me.eq.king.or..not.out1file)
2353 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2354 & ipair(2,i,iset1), qinpair(i,iset1)
2357 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2358 & wfrag_back(3,i,iset1),
2359 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2360 if(me.eq.king.or..not.out1file)
2361 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2362 & wfrag_back(2,i,iset1),
2363 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2364 & ifrag_back(2,i,iset1)
2369 C---------------------------------------------------------------------------
2370 subroutine read_afminp
2371 implicit real*8 (a-h,o-z)
2372 include 'DIMENSIONS'
2376 include 'COMMON.SETUP'
2377 include 'COMMON.CONTROL'
2378 include 'COMMON.CHAIN'
2379 include 'COMMON.IOUNITS'
2380 include 'COMMON.SBRIDGE'
2381 character*320 afmcard
2382 c print *, "wchodze"
2383 call card_concat(afmcard)
2384 call readi(afmcard,"BEG",afmbeg,0)
2385 call readi(afmcard,"END",afmend,0)
2386 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2387 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2388 print *,'FORCE=' ,forceAFMconst
2389 CCCC NOW PROPERTIES FOR AFM
2392 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2394 distafminit=dsqrt(distafminit)
2395 print *,'initdist',distafminit
2399 c-------------------------------------------------------------------------------
2400 subroutine read_dist_constr
2401 implicit real*8 (a-h,o-z)
2402 include 'DIMENSIONS'
2406 include 'COMMON.SETUP'
2407 include 'COMMON.CONTROL'
2408 include 'COMMON.CHAIN'
2409 include 'COMMON.IOUNITS'
2410 include 'COMMON.SBRIDGE'
2411 integer ifrag_(2,100),ipair_(2,100)
2412 double precision wfrag_(100),wpair_(100)
2413 character*500 controlcard
2414 c write (iout,*) "Calling read_dist_constr"
2415 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2417 call card_concat(controlcard)
2418 call readi(controlcard,"NFRAG",nfrag_,0)
2419 call readi(controlcard,"NPAIR",npair_,0)
2420 call readi(controlcard,"NDIST",ndist_,0)
2421 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2422 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2423 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2424 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2425 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2426 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2427 c write (iout,*) "IFRAG"
2429 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2431 c write (iout,*) "IPAIR"
2433 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2435 if (.not.refstr .and. nfrag.gt.0) then
2437 & "ERROR: no reference structure to compute distance restraints"
2439 & "Restraints must be specified explicitly (NDIST=number)"
2442 if (nfrag.lt.2 .and. npair.gt.0) then
2443 write (iout,*) "ERROR: Less than 2 fragments specified",
2444 & " but distance restraints between pairs requested"
2449 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2450 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2451 & ifrag_(2,i)=nstart_sup+nsup-1
2452 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2454 if (wfrag_(i).gt.0.0d0) then
2455 do j=ifrag_(1,i),ifrag_(2,i)-1
2456 do k=j+1,ifrag_(2,i)
2457 c write (iout,*) "j",j," k",k
2459 if (constr_dist.eq.1) then
2464 forcon(nhpb)=wfrag_(i)
2465 else if (constr_dist.eq.2) then
2466 if (ddjk.le.dist_cut) then
2471 forcon(nhpb)=wfrag_(i)
2478 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2481 if (.not.out1file .or. me.eq.king)
2482 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2483 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2485 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2486 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2493 if (wpair_(i).gt.0.0d0) then
2501 do j=ifrag_(1,ii),ifrag_(2,ii)
2502 do k=ifrag_(1,jj),ifrag_(2,jj)
2506 forcon(nhpb)=wpair_(i)
2507 dhpb(nhpb)=dist(j,k)
2509 if (.not.out1file .or. me.eq.king)
2510 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2511 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2513 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2514 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2521 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2522 if (forcon(nhpb+1).gt.0.0d0) then
2524 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2526 if (.not.out1file .or. me.eq.king)
2527 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2528 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2530 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2531 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2538 c-------------------------------------------------------------------------------
2540 subroutine read_constr_homology
2542 include 'DIMENSIONS'
2546 include 'COMMON.SETUP'
2547 include 'COMMON.CONTROL'
2548 include 'COMMON.CHAIN'
2549 include 'COMMON.IOUNITS'
2551 include 'COMMON.GEO'
2552 include 'COMMON.INTERACT'
2553 include 'COMMON.NAMES'
2555 c For new homol impl
2557 include 'COMMON.VAR'
2560 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2562 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2563 c & sigma_odl_temp(maxres,maxres,max_template)
2565 character*24 model_ki_dist, model_ki_angle
2566 character*500 controlcard
2567 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2568 logical lprn /.true./
2573 c FP - Nov. 2014 Temporary specifications for new vars
2575 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2577 double precision, dimension (max_template,maxres) :: rescore
2578 double precision, dimension (max_template,maxres) :: rescore2
2579 double precision, dimension (max_template,maxres) :: rescore3
2580 character*24 pdbfile,tpl_k_rescore
2581 c -----------------------------------------------------------------
2582 c Reading multiple PDB ref structures and calculation of retraints
2583 c not using pre-computed ones stored in files model_ki_{dist,angle}
2585 c -----------------------------------------------------------------
2588 c Alternative: reading from input
2589 call card_concat(controlcard)
2590 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2591 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2592 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2593 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2594 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2595 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2596 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2597 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2598 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2599 if(.not.read2sigma.and.start_from_model) then
2600 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2601 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2602 start_from_model=.false.
2604 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2605 & write(iout,*) 'START_FROM_MODELS is ON'
2606 if(start_from_model .and. rest) then
2607 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2608 write(iout,*) 'START_FROM_MODELS is OFF'
2609 write(iout,*) 'remove restart keyword from input'
2612 if (homol_nset.gt.1)then
2613 call card_concat(controlcard)
2614 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2615 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2616 write(iout,*) "iset homology_weight "
2618 write(iout,*) i,waga_homology(i)
2621 iset=mod(kolor,homol_nset)+1
2624 waga_homology(1)=1.0
2627 cd write (iout,*) "nnt",nnt," nct",nct
2634 c write(iout,*) 'nnt=',nnt,'nct=',nct
2637 do k=1,constr_homology
2650 do k=1,constr_homology
2652 read(inp,'(a)') pdbfile
2653 if(me.eq.king .or. .not. out1file)
2654 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2655 & pdbfile(:ilen(pdbfile))
2656 open(ipdbin,file=pdbfile,status='old',err=33)
2658 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2659 & pdbfile(:ilen(pdbfile))
2662 c print *,'Begin reading pdb data'
2664 c Files containing res sim or local scores (former containing sigmas)
2667 write(kic2,'(bz,i2.2)') k
2669 tpl_k_rescore="template"//kic2//".sco"
2672 if (read2sigma) then
2673 call readpdb_template(k)
2678 c Distance restraints
2681 C Copy the coordinates from reference coordinates (?)
2685 c write (iout,*) "c(",j,i,") =",c(j,i)
2689 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2691 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2692 open (ientin,file=tpl_k_rescore,status='old')
2693 if (nnt.gt.1) rescore(k,1)=0.0d0
2694 do irec=nnt,nct ! loop for reading res sim
2695 if (read2sigma) then
2696 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2697 & rescore3_tmp,idomain_tmp
2699 idomain(k,i_tmp)=idomain_tmp
2700 rescore(k,i_tmp)=rescore_tmp
2701 rescore2(k,i_tmp)=rescore2_tmp
2702 rescore3(k,i_tmp)=rescore3_tmp
2703 if (.not. out1file .or. me.eq.king)
2704 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
2705 & i_tmp,rescore2_tmp,rescore_tmp,
2706 & rescore3_tmp,idomain_tmp
2709 read (ientin,*,end=1401) rescore_tmp
2711 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2712 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2713 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2718 if (waga_dist.ne.0.0d0) then
2726 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2727 c write (iout,*) k,i,j,distal,dist2_cut
2729 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2730 & .and. distal.le.dist2_cut ) then
2736 c write (iout,*) "k",k
2737 c write (iout,*) "i",i," j",j," constr_homology",
2742 if (read2sigma) then
2745 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2747 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2748 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2749 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2751 if (odl(k,ii).le.dist_cut) then
2752 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2755 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2756 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2758 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2759 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2763 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2766 l_homo(k,ii)=.false.
2773 c Theta, dihedral and SC retraints
2775 if (waga_angle.gt.0.0d0) then
2776 c open (ientin,file=tpl_k_sigma_dih,status='old')
2777 c do irec=1,maxres-3 ! loop for reading sigma_dih
2778 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2779 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2780 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2781 c & sigma_dih(k,i+nnt-1)
2786 if (idomain(k,i).eq.0) then
2790 dih(k,i)=phiref(i) ! right?
2791 c read (ientin,*) sigma_dih(k,i) ! original variant
2792 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2793 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2794 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2795 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2796 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2798 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2799 & rescore(k,i-2)+rescore(k,i-3))/4.0
2800 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2801 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2802 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2803 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2804 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2805 c Instead of res sim other local measure of b/b str reliability possible
2806 if (sigma_dih(k,i).ne.0)
2807 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2808 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2813 if (waga_theta.gt.0.0d0) then
2814 c open (ientin,file=tpl_k_sigma_theta,status='old')
2815 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2816 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2817 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2818 c & sigma_theta(k,i+nnt-1)
2823 do i = nnt+2,nct ! right? without parallel.
2824 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2825 c do i=ithet_start,ithet_end ! with FG parallel.
2826 if (idomain(k,i).eq.0) then
2827 sigma_theta(k,i)=0.0
2830 thetatpl(k,i)=thetaref(i)
2831 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2832 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2833 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2834 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2835 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2836 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2837 & rescore(k,i-2))/3.0
2838 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2839 if (sigma_theta(k,i).ne.0)
2840 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2842 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2843 c rescore(k,i-2) ! right expression ?
2844 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2848 if (waga_d.gt.0.0d0) then
2849 c open (ientin,file=tpl_k_sigma_d,status='old')
2850 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2851 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2852 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2853 c & sigma_d(k,i+nnt-1)
2857 do i = nnt,nct ! right? without parallel.
2858 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2859 c do i=loc_start,loc_end ! with FG parallel.
2860 if (itype(i).eq.10) cycle
2861 if (idomain(k,i).eq.0 ) then
2868 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2869 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2870 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2871 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
2872 sigma_d(k,i)=rescore3(k,i) ! right expression ?
2873 if (sigma_d(k,i).ne.0)
2874 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2876 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
2877 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2878 c read (ientin,*) sigma_d(k,i) ! 1st variant
2883 c remove distance restraints not used in any model from the list
2884 c shift data in all arrays
2886 if (waga_dist.ne.0.0d0) then
2892 if (ii_in_use(ii).eq.0.and.liiflag) then
2896 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
2897 & .not.liiflag.and.ii.eq.lim_odl) then
2898 if (ii.eq.lim_odl) then
2899 iishift=ii-iistart+1
2904 do ki=iistart,lim_odl-iishift
2905 ires_homo(ki)=ires_homo(ki+iishift)
2906 jres_homo(ki)=jres_homo(ki+iishift)
2907 ii_in_use(ki)=ii_in_use(ki+iishift)
2908 do k=1,constr_homology
2909 odl(k,ki)=odl(k,ki+iishift)
2910 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
2911 l_homo(k,ki)=l_homo(k,ki+iishift)
2915 lim_odl=lim_odl-iishift
2920 if (constr_homology.gt.0) call homology_partition
2921 if (constr_homology.gt.0) call init_int_table
2922 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2923 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2927 if (.not.lprn) return
2928 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2929 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2930 write (iout,*) "Distance restraints from templates"
2932 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
2933 & ii,ires_homo(ii),jres_homo(ii),
2934 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
2935 & ki=1,constr_homology)
2937 write (iout,*) "Dihedral angle restraints from templates"
2939 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
2940 & (rad2deg*dih(ki,i),
2941 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2943 write (iout,*) "Virtual-bond angle restraints from templates"
2945 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
2946 & (rad2deg*thetatpl(ki,i),
2947 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2949 write (iout,*) "SC restraints from templates"
2951 write(iout,'(i5,100(4f8.2,4x))') i,
2952 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
2953 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2956 c -----------------------------------------------------------------
2959 c----------------------------------------------------------------------
2962 subroutine flush(iu)
2967 subroutine flush(iu)
2972 c------------------------------------------------------------------------------
2973 subroutine copy_to_tmp(source)
2974 include "DIMENSIONS"
2975 include "COMMON.IOUNITS"
2976 character*(*) source
2977 character* 256 tmpfile
2981 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2982 inquire(file=tmpfile,exist=ex)
2984 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2985 & " to temporary directory..."
2986 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2987 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2991 c------------------------------------------------------------------------------
2992 subroutine move_from_tmp(source)
2993 include "DIMENSIONS"
2994 include "COMMON.IOUNITS"
2995 character*(*) source
2998 write (*,*) "Moving ",source(:ilen(source)),
2999 & " from temporary directory to working directory"
3000 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3001 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3004 c------------------------------------------------------------------------------
3005 subroutine random_init(seed)
3007 C Initialize random number generator
3009 implicit real*8 (a-h,o-z)
3010 include 'DIMENSIONS'
3016 logical OKRandom, prng_restart
3018 integer iseed_array(4)
3020 include 'COMMON.IOUNITS'
3021 include 'COMMON.TIME1'
3022 include 'COMMON.THREAD'
3023 include 'COMMON.SBRIDGE'
3024 include 'COMMON.CONTROL'
3025 include 'COMMON.MCM'
3026 include 'COMMON.MAP'
3027 include 'COMMON.HEADER'
3028 include 'COMMON.CSA'
3029 include 'COMMON.CHAIN'
3030 include 'COMMON.MUCA'
3032 include 'COMMON.FFIELD'
3033 include 'COMMON.SETUP'
3034 iseed=-dint(dabs(seed))
3035 if (iseed.eq.0) then
3036 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3037 & 'Random seed undefined. The program will stop.'
3038 write (*,'(/80(1h*)/20x,a/80(1h*))')
3039 & 'Random seed undefined. The program will stop.'
3041 call mpi_finalize(mpi_comm_world,ierr)
3043 stop 'Bad random seed.'
3046 if (fg_rank.eq.0) then
3050 if(me.eq.king .or. .not. out1file)
3051 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3052 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3053 OKRandom = prng_restart(me,iseedi8)
3056 tmp=65536.0d0**(4-i)
3057 iseed_array(i) = dint(seed/tmp)
3058 seed=seed-iseed_array(i)*tmp
3061 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3062 & (iseed_array(i),i=1,4)
3063 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3064 & (iseed_array(i),i=1,4)
3065 OKRandom = prng_restart(me,iseed_array)
3068 c r1 = prng_next(me)
3069 r1=ran_number(0.0D0,1.0D0)
3070 if(me.eq.king .or. .not. out1file)
3071 & write (iout,*) 'ran_num',r1
3072 if (r1.lt.0.0d0) OKRandom=.false.
3074 if (.not.OKRandom) then
3075 write (iout,*) 'PRNG IS NOT WORKING!!!'
3076 print *,'PRNG IS NOT WORKING!!!'
3079 call mpi_abort(mpi_comm_world,error_msg,ierr)
3082 write (iout,*) 'too many processors for parallel prng'
3083 write (*,*) 'too many processors for parallel prng'
3091 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)