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/4a5,2a8,3a10,a5)')
46 & "The following",nhpb-nss,
47 & " distance restraints have been imposed:",
48 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
51 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
52 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
59 c print *,"Processor",myrank," leaves READRTNS"
62 C-------------------------------------------------------------------------------
63 subroutine read_control
67 implicit real*8 (a-h,o-z)
71 logical OKRandom, prng_restart
74 include 'COMMON.IOUNITS'
75 include 'COMMON.TIME1'
76 include 'COMMON.THREAD'
77 include 'COMMON.SBRIDGE'
78 include 'COMMON.CONTROL'
81 include 'COMMON.HEADER'
83 include 'COMMON.CHAIN'
86 include 'COMMON.FFIELD'
87 include 'COMMON.INTERACT'
88 include 'COMMON.SETUP'
89 include 'COMMON.SPLITELE'
90 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
91 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
93 character*320 controlcard
98 read (INP,'(a)') titel
99 call card_concat(controlcard)
100 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
101 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
102 call reada(controlcard,'SEED',seed,0.0D0)
103 call random_init(seed)
104 C Set up the time limit (caution! The time must be input in minutes!)
105 read_cart=index(controlcard,'READ_CART').gt.0
106 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
107 write (iout,*) "constr_dist",constr_dist
108 call readi(controlcard,'NSAXS',nsaxs,0)
109 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
110 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
111 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
112 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
113 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
114 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
115 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
116 call readi(controlcard,'SYM',symetr,1)
117 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
118 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
119 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
120 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
121 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
122 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
123 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
124 call reada(controlcard,'DRMS',drms,0.1D0)
125 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
126 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
127 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
128 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
129 write (iout,'(a,f10.1)')'DRMS = ',drms
130 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
131 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
133 call readi(controlcard,'NZ_START',nz_start,0)
134 call readi(controlcard,'NZ_END',nz_end,0)
135 call readi(controlcard,'IZ_SC',iz_sc,0)
137 safety = 60.0d0*safety
140 call reada(controlcard,"T_BATH",t_bath,300.0d0)
141 minim=(index(controlcard,'MINIMIZE').gt.0)
142 dccart=(index(controlcard,'CART').gt.0)
143 overlapsc=(index(controlcard,'OVERLAP').gt.0)
144 overlapsc=.not.overlapsc
145 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
146 searchsc=.not.searchsc
147 sideadd=(index(controlcard,'SIDEADD').gt.0)
148 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
149 outpdb=(index(controlcard,'PDBOUT').gt.0)
150 outmol2=(index(controlcard,'MOL2OUT').gt.0)
151 pdbref=(index(controlcard,'PDBREF').gt.0)
152 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
153 indpdb=index(controlcard,'PDBSTART')
154 extconf=(index(controlcard,'EXTCONF').gt.0)
155 AFMlog=(index(controlcard,'AFM'))
156 selfguide=(index(controlcard,'SELFGUIDE'))
157 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
158 call readi(controlcard,'IPRINT',iprint,0)
159 call readi(controlcard,'MAXGEN',maxgen,10000)
160 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
161 call readi(controlcard,"KDIAG",kdiag,0)
162 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
163 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
164 & write (iout,*) "RESCALE_MODE",rescale_mode
165 split_ene=index(controlcard,'SPLIT_ENE').gt.0
166 if (index(controlcard,'REGULAR').gt.0.0D0) then
167 call reada(controlcard,'WEIDIS',weidis,0.1D0)
171 if (index(controlcard,'CHECKGRAD').gt.0) then
173 if (index(controlcard,'CART').gt.0) then
175 elseif (index(controlcard,'CARINT').gt.0) then
180 elseif (index(controlcard,'THREAD').gt.0) then
182 call readi(controlcard,'THREAD',nthread,0)
183 if (nthread.gt.0) then
184 call reada(controlcard,'WEIDIS',weidis,0.1D0)
187 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
188 stop 'Error termination in Read_Control.'
190 else if (index(controlcard,'MCMA').gt.0) then
192 else if (index(controlcard,'MCEE').gt.0) then
194 else if (index(controlcard,'MULTCONF').gt.0) then
196 else if (index(controlcard,'MAP').gt.0) then
198 call readi(controlcard,'MAP',nmap,0)
199 else if (index(controlcard,'CSA').gt.0) then
201 crc else if (index(controlcard,'ZSCORE').gt.0) then
203 crc ZSCORE is rm from UNRES, modecalc=9 is available
206 cfcm else if (index(controlcard,'MCMF').gt.0) then
208 else if (index(controlcard,'SOFTREG').gt.0) then
210 else if (index(controlcard,'CHECK_BOND').gt.0) then
212 else if (index(controlcard,'TEST').gt.0) then
214 else if (index(controlcard,'MD').gt.0) then
216 else if (index(controlcard,'RE ').gt.0) then
220 lmuca=index(controlcard,'MUCA').gt.0
221 call readi(controlcard,'MUCADYN',mucadyn,0)
222 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
223 if (lmuca .and. (me.eq.king .or. .not.out1file ))
225 write (iout,*) 'MUCADYN=',mucadyn
226 write (iout,*) 'MUCASMOOTH=',muca_smooth
229 iscode=index(controlcard,'ONE_LETTER')
230 indphi=index(controlcard,'PHI')
231 indback=index(controlcard,'BACK')
232 iranconf=index(controlcard,'RAND_CONF')
233 i2ndstr=index(controlcard,'USE_SEC_PRED')
234 gradout=index(controlcard,'GRADOUT').gt.0
235 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
236 C DISTCHAINMAX become obsolete for periodic boundry condition
237 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
238 C Reading the dimensions of box in x,y,z coordinates
239 call reada(controlcard,'BOXX',boxxsize,100.0d0)
240 call reada(controlcard,'BOXY',boxysize,100.0d0)
241 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
242 c Cutoff range for interactions
243 call reada(controlcard,"R_CUT",r_cut,15.0d0)
244 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
245 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
246 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
247 if (lipthick.gt.0.0d0) then
248 bordliptop=(boxzsize+lipthick)/2.0
249 bordlipbot=bordliptop-lipthick
251 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
252 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
253 buflipbot=bordlipbot+lipbufthick
254 bufliptop=bordliptop-lipbufthick
255 if ((lipbufthick*2.0d0).gt.lipthick)
256 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
258 c write(iout,*) "bordliptop=",bordliptop
259 c write(iout,*) "bordlipbot=",bordlipbot
260 c write(iout,*) "bufliptop=",bufliptop
261 c write(iout,*) "buflipbot=",buflipbot
264 if (me.eq.king .or. .not.out1file )
265 & write (iout,*) "DISTCHAINMAX",distchainmax
267 if(me.eq.king.or..not.out1file)
268 & write (iout,'(2a)') diagmeth(kdiag),
269 & ' routine used to diagonalize matrices.'
272 c--------------------------------------------------------------------------
273 subroutine read_REMDpar
277 implicit real*8 (a-h,o-z)
279 include 'COMMON.IOUNITS'
280 include 'COMMON.TIME1'
283 include 'COMMON.LANGEVIN'
285 include 'COMMON.LANGEVIN.lang0'
287 include 'COMMON.INTERACT'
288 include 'COMMON.NAMES'
290 include 'COMMON.REMD'
291 include 'COMMON.CONTROL'
292 include 'COMMON.SETUP'
294 character*320 controlcard
295 character*3200 controlcard1
296 integer iremd_m_total
298 if(me.eq.king.or..not.out1file)
299 & write (iout,*) "REMD setup"
301 call card_concat(controlcard)
302 call readi(controlcard,"NREP",nrep,3)
303 call readi(controlcard,"NSTEX",nstex,1000)
304 call reada(controlcard,"RETMIN",retmin,10.0d0)
305 call reada(controlcard,"RETMAX",retmax,1000.0d0)
306 mremdsync=(index(controlcard,'SYNC').gt.0)
307 call readi(controlcard,"NSYN",i_sync_step,100)
308 restart1file=(index(controlcard,'REST1FILE').gt.0)
309 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
310 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
311 if(max_cache_traj_use.gt.max_cache_traj)
312 & max_cache_traj_use=max_cache_traj
313 if(me.eq.king.or..not.out1file) then
314 cd if (traj1file) then
315 crc caching is in testing - NTWX is not ignored
316 cd write (iout,*) "NTWX value is ignored"
317 cd write (iout,*) " trajectory is stored to one file by master"
318 cd write (iout,*) " before exchange at NSTEX intervals"
320 write (iout,*) "NREP= ",nrep
321 write (iout,*) "NSTEX= ",nstex
322 write (iout,*) "SYNC= ",mremdsync
323 write (iout,*) "NSYN= ",i_sync_step
324 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
327 if (index(controlcard,'TLIST').gt.0) then
329 call card_concat(controlcard1)
330 read(controlcard1,*) (remd_t(i),i=1,nrep)
331 if(me.eq.king.or..not.out1file)
332 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
335 if (index(controlcard,'MLIST').gt.0) then
337 call card_concat(controlcard1)
338 read(controlcard1,*) (remd_m(i),i=1,nrep)
339 if(me.eq.king.or..not.out1file) then
340 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
343 iremd_m_total=iremd_m_total+remd_m(i)
345 write (iout,*) 'Total number of replicas ',iremd_m_total
348 if(me.eq.king.or..not.out1file)
349 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
352 c--------------------------------------------------------------------------
353 subroutine read_MDpar
357 implicit real*8 (a-h,o-z)
359 include 'COMMON.IOUNITS'
360 include 'COMMON.TIME1'
363 include 'COMMON.LANGEVIN'
365 include 'COMMON.LANGEVIN.lang0'
367 include 'COMMON.INTERACT'
368 include 'COMMON.NAMES'
370 include 'COMMON.SETUP'
371 include 'COMMON.CONTROL'
372 include 'COMMON.SPLITELE'
374 character*320 controlcard
376 call card_concat(controlcard)
377 call readi(controlcard,"NSTEP",n_timestep,1000000)
378 call readi(controlcard,"NTWE",ntwe,100)
379 call readi(controlcard,"NTWX",ntwx,1000)
380 call reada(controlcard,"DT",d_time,1.0d-1)
381 call reada(controlcard,"DVMAX",dvmax,2.0d1)
382 call reada(controlcard,"DAMAX",damax,1.0d1)
383 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
384 call readi(controlcard,"LANG",lang,0)
385 RESPA = index(controlcard,"RESPA") .gt. 0
386 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
387 ntime_split0=ntime_split
388 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
389 ntime_split0=ntime_split
390 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
391 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
392 rest = index(controlcard,"REST").gt.0
393 tbf = index(controlcard,"TBF").gt.0
394 usampl = index(controlcard,"USAMPL").gt.0
396 mdpdb = index(controlcard,"MDPDB").gt.0
397 call reada(controlcard,"T_BATH",t_bath,300.0d0)
398 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
399 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
400 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
401 if (count_reset_moment.eq.0) count_reset_moment=1000000000
402 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
403 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
404 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
405 if (count_reset_vel.eq.0) count_reset_vel=1000000000
406 large = index(controlcard,"LARGE").gt.0
407 print_compon = index(controlcard,"PRINT_COMPON").gt.0
408 rattle = index(controlcard,"RATTLE").gt.0
409 preminim = index(controlcard,"PREMINIM").gt.0
411 dccart=(index(controlcard,'CART').gt.0)
414 c if performing umbrella sampling, fragments constrained are read from the fragment file
420 if(me.eq.king.or..not.out1file) then
422 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
424 write (iout,'(a)') "The units are:"
425 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
426 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
427 & " acceleration: angstrom/(48.9 fs)**2"
428 write (iout,'(a)') "energy: kcal/mol, temperature: K"
430 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
431 write (iout,'(a60,f10.5,a)')
432 & "Initial time step of numerical integration:",d_time,
434 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
436 write (iout,'(2a,i4,a)')
437 & "A-MTS algorithm used; initial time step for fast-varying",
438 & " short-range forces split into",ntime_split," steps."
439 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
440 & r_cut," lambda",rlamb
442 write (iout,'(2a,f10.5)')
443 & "Maximum acceleration threshold to reduce the time step",
444 & "/increase split number:",damax
445 write (iout,'(2a,f10.5)')
446 & "Maximum predicted energy drift to reduce the timestep",
447 & "/increase split number:",edriftmax
448 write (iout,'(a60,f10.5)')
449 & "Maximum velocity threshold to reduce velocities:",dvmax
450 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
451 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
452 if (rattle) write (iout,'(a60)')
453 & "Rattle algorithm used to constrain the virtual bonds"
454 if (preminim .or. iranconf.gt.0) then
456 & "Initial structure will be energy-minimized"
461 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
462 call reada(controlcard,"RWAT",rwat,1.4d0)
463 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
464 surfarea=index(controlcard,"SURFAREA").gt.0
465 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
466 if(me.eq.king.or..not.out1file)then
467 write (iout,'(/a,$)') "Langevin dynamics calculation"
470 & " with direct integration of Langevin equations"
471 else if (lang.eq.2) then
472 write (iout,'(a/)') " with TINKER stochasic MD integrator"
473 else if (lang.eq.3) then
474 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
475 else if (lang.eq.4) then
476 write (iout,'(a/)') " in overdamped mode"
478 write (iout,'(//a,i5)')
479 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
482 write (iout,'(a60,f10.5)') "Temperature:",t_bath
483 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
484 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
485 write (iout,'(a60,f10.5)')
486 & "Scaling factor of the friction forces:",scal_fric
487 if (surfarea) write (iout,'(2a,i10,a)')
488 & "Friction coefficients will be scaled by solvent-accessible",
489 & " surface area every",reset_fricmat," steps."
491 c Calculate friction coefficients and bounds of stochastic forces
492 eta=6*pi*cPoise*etawat
493 if(me.eq.king.or..not.out1file)
494 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
496 gamp=scal_fric*(pstok+rwat)*eta
497 stdfp=dsqrt(2*Rb*t_bath/d_time)
499 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
500 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
502 if(me.eq.king.or..not.out1file)then
503 write (iout,'(/2a/)')
504 & "Radii of site types and friction coefficients and std's of",
505 & " stochastic forces of fully exposed sites"
506 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
508 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
509 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
513 if(me.eq.king.or..not.out1file)then
514 write (iout,'(a)') "Berendsen bath calculation"
515 write (iout,'(a60,f10.5)') "Temperature:",t_bath
516 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
518 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
519 & count_reset_moment," steps"
521 & write (iout,'(a,i10,a)')
522 & "Velocities will be reset at random every",count_reset_vel,
526 if(me.eq.king.or..not.out1file)
527 & write (iout,'(a31)') "Microcanonical mode calculation"
529 if(me.eq.king.or..not.out1file)then
530 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
532 write(iout,*) "MD running with constraints."
533 write(iout,*) "Equilibration time ", eq_time, " mtus."
534 write(iout,*) "Constraining ", nfrag," fragments."
535 write(iout,*) "Length of each fragment, weight and q0:"
537 write (iout,*) "Set of restraints #",iset
539 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
540 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
542 write(iout,*) "constraints between ", npair, "fragments."
543 write(iout,*) "constraint pairs, weights and q0:"
545 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
546 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
548 write(iout,*) "angle constraints within ", nfrag_back,
549 & "backbone fragments."
550 write(iout,*) "fragment, weights:"
552 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
553 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
554 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
557 iset=mod(kolor,nset)+1
560 if(me.eq.king.or..not.out1file)
561 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
564 c------------------------------------------------------------------------------
567 C Read molecular data.
569 implicit real*8 (a-h,o-z)
575 include 'COMMON.IOUNITS'
578 include 'COMMON.INTERACT'
579 include 'COMMON.LOCAL'
580 include 'COMMON.NAMES'
581 include 'COMMON.CHAIN'
582 include 'COMMON.FFIELD'
583 include 'COMMON.SBRIDGE'
584 include 'COMMON.HEADER'
585 include 'COMMON.CONTROL'
586 include 'COMMON.DBASE'
587 include 'COMMON.THREAD'
588 include 'COMMON.CONTACTS'
589 include 'COMMON.TORCNSTR'
590 include 'COMMON.TIME1'
591 include 'COMMON.BOUNDS'
593 include 'COMMON.SETUP'
594 character*4 sequence(maxres)
596 double precision x(maxvar)
597 character*256 pdbfile
598 character*320 weightcard
599 character*80 weightcard_t,ucase
600 dimension itype_pdb(maxres)
601 common /pizda/ itype_pdb
602 logical seq_comp,fail
603 double precision energia(0:n_ene)
609 C Read weights of the subsequent energy terms.
610 call card_concat(weightcard)
611 call reada(weightcard,'WLONG',wlong,1.0D0)
612 call reada(weightcard,'WSC',wsc,wlong)
613 call reada(weightcard,'WSCP',wscp,wlong)
614 call reada(weightcard,'WELEC',welec,1.0D0)
615 call reada(weightcard,'WVDWPP',wvdwpp,welec)
616 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
617 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
618 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
619 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
620 call reada(weightcard,'WTURN3',wturn3,1.0D0)
621 call reada(weightcard,'WTURN4',wturn4,1.0D0)
622 call reada(weightcard,'WTURN6',wturn6,1.0D0)
623 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
624 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
625 call reada(weightcard,'WBOND',wbond,1.0D0)
626 call reada(weightcard,'WTOR',wtor,1.0D0)
627 call reada(weightcard,'WTORD',wtor_d,1.0D0)
628 call reada(weightcard,'WANG',wang,1.0D0)
629 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
630 call reada(weightcard,'SCAL14',scal14,0.4D0)
631 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
632 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
633 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
634 call reada(weightcard,'TEMP0',temp0,300.0d0)
635 call reada(weightcard,'WLT',wliptran,0.0D0)
636 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
637 if (index(weightcard,'SOFT').gt.0) ipot=6
638 C 12/1/95 Added weight for the multi-body term WCORR
639 call reada(weightcard,'WCORRH',wcorr,1.0D0)
640 if (wcorr4.gt.0.0d0) wcorr=wcorr4
661 if(me.eq.king.or..not.out1file)
662 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
663 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
665 10 format (/'Energy-term weights (unscaled):'//
666 & 'WSCC= ',f10.6,' (SC-SC)'/
667 & 'WSCP= ',f10.6,' (SC-p)'/
668 & 'WELEC= ',f10.6,' (p-p electr)'/
669 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
670 & 'WBOND= ',f10.6,' (stretching)'/
671 & 'WANG= ',f10.6,' (bending)'/
672 & 'WSCLOC= ',f10.6,' (SC local)'/
673 & 'WTOR= ',f10.6,' (torsional)'/
674 & 'WTORD= ',f10.6,' (double torsional)'/
675 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
676 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
677 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
678 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
679 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
680 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
681 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
682 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
683 & 'WTURN6= ',f10.6,' (turns, 6th order)')
684 if(me.eq.king.or..not.out1file)then
685 if (wcorr4.gt.0.0d0) then
686 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
687 & 'between contact pairs of peptide groups'
688 write (iout,'(2(a,f5.3/))')
689 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
690 & 'Range of quenching the correlation terms:',2*delt_corr
691 else if (wcorr.gt.0.0d0) then
692 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
693 & 'between contact pairs of peptide groups'
695 write (iout,'(a,f8.3)')
696 & 'Scaling factor of 1,4 SC-p interactions:',scal14
697 write (iout,'(a,f8.3)')
698 & 'General scaling factor of SC-p interactions:',scalscp
700 r0_corr=cutoff_corr-delt_corr
702 aad(i,1)=scalscp*aad(i,1)
703 aad(i,2)=scalscp*aad(i,2)
704 bad(i,1)=scalscp*bad(i,1)
705 bad(i,2)=scalscp*bad(i,2)
707 call rescale_weights(t_bath)
708 if(me.eq.king.or..not.out1file)
709 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
710 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
712 22 format (/'Energy-term weights (scaled):'//
713 & 'WSCC= ',f10.6,' (SC-SC)'/
714 & 'WSCP= ',f10.6,' (SC-p)'/
715 & 'WELEC= ',f10.6,' (p-p electr)'/
716 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
717 & 'WBOND= ',f10.6,' (stretching)'/
718 & 'WANG= ',f10.6,' (bending)'/
719 & 'WSCLOC= ',f10.6,' (SC local)'/
720 & 'WTOR= ',f10.6,' (torsional)'/
721 & 'WTORD= ',f10.6,' (double torsional)'/
722 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
723 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
724 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
725 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
726 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
727 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
728 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
729 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
730 & 'WTURN6= ',f10.6,' (turns, 6th order)')
731 if(me.eq.king.or..not.out1file)
732 & write (iout,*) "Reference temperature for weights calculation:",
734 call reada(weightcard,"D0CM",d0cm,3.78d0)
735 call reada(weightcard,"AKCM",akcm,15.1d0)
736 call reada(weightcard,"AKTH",akth,11.0d0)
737 call reada(weightcard,"AKCT",akct,12.0d0)
738 call reada(weightcard,"V1SS",v1ss,-1.08d0)
739 call reada(weightcard,"V2SS",v2ss,7.61d0)
740 call reada(weightcard,"V3SS",v3ss,13.7d0)
741 call reada(weightcard,"EBR",ebr,-5.50D0)
742 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
744 dyn_ss_mask(i)=.false.
748 dyn_ssbond_ij(i,j)=1.0d300
751 call reada(weightcard,"HT",Ht,0.0D0)
753 ss_depth=ebr/wsc-0.25*eps(1,1)
754 Ht=Ht/wsc-0.25*eps(1,1)
755 akcm=akcm*wstrain/wsc
756 akth=akth*wstrain/wsc
757 akct=akct*wstrain/wsc
758 v1ss=v1ss*wstrain/wsc
759 v2ss=v2ss*wstrain/wsc
760 v3ss=v3ss*wstrain/wsc
762 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
765 if(me.eq.king.or..not.out1file) then
766 write (iout,*) "Parameters of the SS-bond potential:"
767 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
769 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
770 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
771 write (iout,*)" HT",Ht
772 print *,'indpdb=',indpdb,' pdbref=',pdbref
774 if (indpdb.gt.0 .or. pdbref) then
775 read(inp,'(a)') pdbfile
776 if(me.eq.king.or..not.out1file)
777 & write (iout,'(2a)') 'PDB data will be read from file ',
778 & pdbfile(:ilen(pdbfile))
779 open(ipdbin,file=pdbfile,status='old',err=33)
781 33 write (iout,'(a)') 'Error opening PDB file.'
784 c print *,'Begin reading pdb data'
793 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
794 & (crefjlee(j,i+nres),j=1,3)
797 c print *,'Finished reading pdb data'
798 if(me.eq.king.or..not.out1file)
799 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
800 & ' nstart_sup=',nstart_sup
802 itype_pdb(i)=itype(i)
806 nct=nstart_sup+nsup-1
807 call contact(.false.,ncont_ref,icont_ref,co)
810 if(me.eq.king.or..not.out1file)
811 & write(iout,*)'Adding sidechains'
815 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
818 do while (fail.and.nsi.le.maxsi)
819 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
823 c Calculalte only the coordinates of the current sidechain; no need to rebuild
825 call locate_side_chain(i)
826 if(fail) write(iout,*)'Adding sidechain failed for res ',
827 & i,' after ',nsi,' trials'
830 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
834 if (indpdb.eq.0) then
835 C Read sequence if not taken from the pdb file.
837 c print *,'nres=',nres
838 if (iscode.gt.0) then
839 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
841 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
843 C Convert sequence to numeric code
845 itype(i)=rescode(i,sequence(i),iscode)
847 C Assign initial virtual bond lengths
853 vbld(i+nres)=dsc(iabs(itype(i)))
854 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
855 c write (iout,*) "i",i," itype",itype(i),
856 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
860 c print '(20i4)',(itype(i),i=1,nres)
863 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
865 if (itype(i).eq.ntyp1) then
869 else if (iabs(itype(i+1)).ne.20) then
871 else if (iabs(itype(i)).ne.20) then
878 if(me.eq.king.or..not.out1file)then
879 write (iout,*) "ITEL"
881 write (iout,*) i,itype(i),itel(i)
883 print *,'Call Read_Bridge.'
886 C 8/13/98 Set limits to generating the dihedral angles
891 read (inp,*) ndih_constr
892 write (iout,*) "ndish_constr",ndih_constr
893 if (ndih_constr.gt.0) then
895 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
896 if(me.eq.king.or..not.out1file)then
898 & 'There are',ndih_constr,' constraints on phi angles.'
900 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
904 phi0(i)=deg2rad*phi0(i)
905 drange(i)=deg2rad*drange(i)
907 if(me.eq.king.or..not.out1file)
908 & write (iout,*) 'FTORS',ftors
911 phibound(1,ii) = phi0(i)-drange(i)
912 phibound(2,ii) = phi0(i)+drange(i)
919 write (iout,'(a)') 'Boundaries in phi angle sampling:'
921 write (iout,'(a3,i5,2f10.1)')
922 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
928 cd print *,'NNT=',NNT,' NCT=',NCT
929 if (itype(1).eq.ntyp1) nnt=2
930 if (itype(nres).eq.ntyp1) nct=nct-1
932 if(me.eq.king.or..not.out1file)
933 & write (iout,'(a,i3)') 'nsup=',nsup
935 if (nsup.le.(nct-nnt+1)) then
936 do i=0,nct-nnt+1-nsup
937 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
943 & 'Error - sequences to be superposed do not match.'
946 do i=0,nsup-(nct-nnt+1)
947 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
949 nstart_sup=nstart_sup+i
955 & 'Error - sequences to be superposed do not match.'
958 if (nsup.eq.0) nsup=nct-nnt
959 if (nstart_sup.eq.0) nstart_sup=nnt
960 if (nstart_seq.eq.0) nstart_seq=nnt
961 if(me.eq.king.or..not.out1file)
962 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
963 & ' nstart_seq=',nstart_seq
965 c--- Zscore rms -------
966 if (nz_start.eq.0) nz_start=nnt
967 if (nz_end.eq.0 .and. nsup.gt.0) then
969 else if (nz_end.eq.0) then
972 if(me.eq.king.or..not.out1file)then
973 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
974 write (iout,*) 'IZ_SC=',iz_sc
976 c----------------------
979 if (.not.pdbref) then
980 call read_angles(inp,*38)
982 38 write (iout,'(a)') 'Error reading reference structure.'
984 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
985 stop 'Error reading reference structure'
987 39 call chainbuild_extconf
989 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
999 call contact(.true.,ncont_ref,icont_ref,co)
1001 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1003 if (constr_dist.gt.0) call read_dist_constr
1004 c write (iout,*) "After read_dist_constr nhpb",nhpb
1005 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1006 if(me.eq.king.or..not.out1file)
1007 & write (iout,*) 'Contact order:',co
1009 if(me.eq.king.or..not.out1file)
1010 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1013 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1015 if(me.eq.king.or..not.out1file)
1016 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1017 & icont_ref(1,i),' ',
1018 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1022 write (iout,*) "calling read_saxs_consrtr",nsaxs
1023 if (nsaxs.gt.0) call read_saxs_constr
1025 if (constr_homology.gt.0) then
1026 call read_constr_homology
1027 if (indpdb.gt.0 .or. pdbref) then
1030 c(j,i)=crefjlee(j,i)
1031 cref(j,i,1)=crefjlee(j,i)
1036 write (iout,*) "Array C"
1038 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1039 & (c(j,i+nres),j=1,3)
1041 write (iout,*) "Array Cref"
1043 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1044 & (cref(j,i+nres,1),j=1,3)
1047 call int_from_cart1(.false.)
1048 call sc_loc_geom(.false.)
1050 thetaref(i)=theta(i)
1055 dc(j,i)=c(j,i+1)-c(j,i)
1056 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1061 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1062 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1070 if (nhpb.gt.0) call hpb_partition
1071 c write (iout,*) "After read_dist_constr nhpb",nhpb
1073 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1074 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1075 & modecalc.ne.10) then
1076 C If input structure hasn't been supplied from the PDB file read or generate
1078 if (iranconf.eq.0 .and. .not. extconf) then
1079 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1080 & write (iout,'(a)') 'Initial geometry will be read in.'
1082 read(inp,'(8f10.5)',end=36,err=36)
1083 & ((c(l,k),l=1,3),k=1,nres),
1084 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1085 write (iout,*) "Exit READ_CART"
1086 write (iout,'(8f10.5)')
1087 & ((c(l,k),l=1,3),k=1,nres),
1088 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1089 call int_from_cart1(.true.)
1090 write (iout,*) "Finish INT_TO_CART"
1093 dc(j,i)=c(j,i+1)-c(j,i)
1094 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1098 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1100 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1101 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1107 write (iout,*) "Calling read_ang"
1108 call read_angles(inp,*36)
1109 write (iout,*) "Calling chainbuild"
1110 call chainbuild_extconf
1113 36 write (iout,'(a)') 'Error reading angle file.'
1115 call mpi_finalize( MPI_COMM_WORLD,IERR )
1117 stop 'Error reading angle file.'
1119 else if (extconf) then
1120 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1121 & write (iout,'(a)') 'Extended chain initial geometry.'
1123 theta(i)=90d0*deg2rad
1126 phi(i)=180d0*deg2rad
1129 alph(i)=110d0*deg2rad
1132 omeg(i)=-120d0*deg2rad
1133 if (itype(i).le.0) omeg(i)=-omeg(i)
1135 c from old chainbuild
1137 C Define the origin and orientation of the coordinate system and locate the
1138 C first three CA's and SC(2).
1142 * Build the alpha-carbon chain.
1145 call locate_next_res(i)
1148 C First and last SC must coincide with the corresponding CA.
1152 dc_norm(j,nres+1)=0.0D0
1153 dc(j,nres+nres)=0.0D0
1154 dc_norm(j,nres+nres)=0.0D0
1156 c(j,nres+nres)=c(j,nres)
1159 C Define the origin and orientation of the coordinate system and locate the
1160 C first three CA's and SC(2).
1164 * Build the alpha-carbon chain.
1167 call locate_next_res(i)
1170 C First and last SC must coincide with the corresponding CA.
1174 dc_norm(j,nres+1)=0.0D0
1175 dc(j,nres+nres)=0.0D0
1176 dc_norm(j,nres+nres)=0.0D0
1178 c(j,nres+nres)=c(j,nres)
1183 if(me.eq.king.or..not.out1file)
1184 & write (iout,'(a)') 'Random-generated initial geometry.'
1188 if (me.eq.king .or. fg_rank.eq.0 .and. (
1189 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1193 call gen_rand_conf(itmp,*30)
1195 30 write (iout,*) 'Failed to generate random conformation',
1196 & ', itrial=',itrial
1197 write (*,*) 'Processor:',me,
1198 & ' Failed to generate random conformation',
1208 write (iout,'(a,i3,a)') 'Processor:',me,
1209 & ' error in generating random conformation.'
1210 write (*,'(a,i3,a)') 'Processor:',me,
1211 & ' error in generating random conformation.'
1214 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1219 & ' error in generating random conformation.'
1224 elseif (modecalc.eq.4) then
1225 read (inp,'(a)') intinname
1226 open (intin,file=intinname,status='old',err=333)
1227 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1228 & write (iout,'(a)') 'intinname',intinname
1229 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1231 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1233 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1235 stop 'Error opening angle file.'
1239 C Generate distance constraints, if the PDB structure is to be regularized.
1240 if (nthread.gt.0) then
1241 call read_threadbase
1244 if (me.eq.king .or. .not. out1file)
1246 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1247 write (iout,'(/a,i3,a)')
1248 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1249 write (iout,'(20i4)') (iss(i),i=1,ns)
1251 write(iout,*)"Running with dynamic disulfide-bond formation"
1253 write (iout,'(/a/)') 'Pre-formed links are:'
1259 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1260 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1266 if (ns.gt.0.and.dyn_ss) then
1270 forcon(i-nss)=forcon(i)
1277 dyn_ss_mask(iss(i))=.true.
1280 if (i2ndstr.gt.0) call secstrp2dihc
1281 c call geom_to_var(nvar,x)
1282 c call etotal(energia(0))
1283 c call enerprint(energia(0))
1284 c call briefout(0,etot)
1286 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1287 cd write (iout,'(a)') 'Variable list:'
1288 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1290 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1291 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1292 & 'Processor',myrank,': end reading molecular data.'
1296 c--------------------------------------------------------------------------
1297 logical function seq_comp(itypea,itypeb,length)
1299 integer length,itypea(length),itypeb(length)
1302 if (itypea(i).ne.itypeb(i)) then
1310 c-----------------------------------------------------------------------------
1311 subroutine read_bridge
1312 C Read information about disulfide bridges.
1313 implicit real*8 (a-h,o-z)
1314 include 'DIMENSIONS'
1318 include 'COMMON.IOUNITS'
1319 include 'COMMON.GEO'
1320 include 'COMMON.VAR'
1321 include 'COMMON.INTERACT'
1322 include 'COMMON.LOCAL'
1323 include 'COMMON.NAMES'
1324 include 'COMMON.CHAIN'
1325 include 'COMMON.FFIELD'
1326 include 'COMMON.SBRIDGE'
1327 include 'COMMON.HEADER'
1328 include 'COMMON.CONTROL'
1329 include 'COMMON.DBASE'
1330 include 'COMMON.THREAD'
1331 include 'COMMON.TIME1'
1332 include 'COMMON.SETUP'
1333 C Read bridging residues.
1334 read (inp,*) ns,(iss(i),i=1,ns)
1336 if(me.eq.king.or..not.out1file)
1337 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1338 C Check whether the specified bridging residues are cystines.
1340 if (itype(iss(i)).ne.1) then
1341 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1342 & 'Do you REALLY think that the residue ',
1343 & restyp(itype(iss(i))),i,
1344 & ' can form a disulfide bridge?!!!'
1345 write (*,'(2a,i3,a)')
1346 & 'Do you REALLY think that the residue ',
1347 & restyp(itype(iss(i))),i,
1348 & ' can form a disulfide bridge?!!!'
1350 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1355 C Read preformed bridges.
1357 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1359 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1362 C Check if the residues involved in bridges are in the specified list of
1363 C bridging residues.
1366 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1367 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1368 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1369 & ' contains residues present in other pairs.'
1370 write (*,'(a,i3,a)') 'Disulfide pair',i,
1371 & ' contains residues present in other pairs.'
1373 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1379 if (ihpb(i).eq.iss(j)) goto 10
1381 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1384 if (jhpb(i).eq.iss(j)) goto 20
1386 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1392 ihpb(i)=ihpb(i)+nres
1393 jhpb(i)=jhpb(i)+nres
1399 c----------------------------------------------------------------------------
1400 subroutine read_x(kanal,*)
1401 implicit real*8 (a-h,o-z)
1402 include 'DIMENSIONS'
1403 include 'COMMON.GEO'
1404 include 'COMMON.VAR'
1405 include 'COMMON.CHAIN'
1406 include 'COMMON.IOUNITS'
1407 include 'COMMON.CONTROL'
1408 include 'COMMON.LOCAL'
1409 include 'COMMON.INTERACT'
1410 c Read coordinates from input
1412 read(kanal,'(8f10.5)',end=10,err=10)
1413 & ((c(l,k),l=1,3),k=1,nres),
1414 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1417 c(j,2*nres)=c(j,nres)
1419 call int_from_cart1(.false.)
1422 dc(j,i)=c(j,i+1)-c(j,i)
1423 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1427 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1429 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1430 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1438 c----------------------------------------------------------------------------
1439 subroutine read_threadbase
1440 implicit real*8 (a-h,o-z)
1441 include 'DIMENSIONS'
1442 include 'COMMON.IOUNITS'
1443 include 'COMMON.GEO'
1444 include 'COMMON.VAR'
1445 include 'COMMON.INTERACT'
1446 include 'COMMON.LOCAL'
1447 include 'COMMON.NAMES'
1448 include 'COMMON.CHAIN'
1449 include 'COMMON.FFIELD'
1450 include 'COMMON.SBRIDGE'
1451 include 'COMMON.HEADER'
1452 include 'COMMON.CONTROL'
1453 include 'COMMON.DBASE'
1454 include 'COMMON.THREAD'
1455 include 'COMMON.TIME1'
1456 C Read pattern database for threading.
1457 read (icbase,*) nseq
1459 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1460 & nres_base(2,i),nres_base(3,i)
1461 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1463 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1464 c & nres_base(2,i),nres_base(3,i)
1465 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1469 if (weidis.eq.0.0D0) weidis=0.1D0
1478 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1479 write (iout,'(a,i5)') 'nexcl: ',nexcl
1480 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1483 c------------------------------------------------------------------------------
1484 subroutine setup_var
1485 implicit real*8 (a-h,o-z)
1486 include 'DIMENSIONS'
1487 include 'COMMON.IOUNITS'
1488 include 'COMMON.GEO'
1489 include 'COMMON.VAR'
1490 include 'COMMON.INTERACT'
1491 include 'COMMON.LOCAL'
1492 include 'COMMON.NAMES'
1493 include 'COMMON.CHAIN'
1494 include 'COMMON.FFIELD'
1495 include 'COMMON.SBRIDGE'
1496 include 'COMMON.HEADER'
1497 include 'COMMON.CONTROL'
1498 include 'COMMON.DBASE'
1499 include 'COMMON.THREAD'
1500 include 'COMMON.TIME1'
1501 C Set up variable list.
1507 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1509 ialph(i,1)=nvar+nside
1513 if (indphi.gt.0) then
1515 else if (indback.gt.0) then
1520 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1523 c----------------------------------------------------------------------------
1524 subroutine gen_dist_constr
1525 C Generate CA distance constraints.
1526 implicit real*8 (a-h,o-z)
1527 include 'DIMENSIONS'
1528 include 'COMMON.IOUNITS'
1529 include 'COMMON.GEO'
1530 include 'COMMON.VAR'
1531 include 'COMMON.INTERACT'
1532 include 'COMMON.LOCAL'
1533 include 'COMMON.NAMES'
1534 include 'COMMON.CHAIN'
1535 include 'COMMON.FFIELD'
1536 include 'COMMON.SBRIDGE'
1537 include 'COMMON.HEADER'
1538 include 'COMMON.CONTROL'
1539 include 'COMMON.DBASE'
1540 include 'COMMON.THREAD'
1541 include 'COMMON.TIME1'
1542 dimension itype_pdb(maxres)
1543 common /pizda/ itype_pdb
1545 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1546 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1547 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1549 do i=nstart_sup,nstart_sup+nsup-1
1550 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1551 cd & ' seq_pdb', restyp(itype_pdb(i))
1552 do j=i+2,nstart_sup+nsup-1
1554 ihpb(nhpb)=i+nstart_seq-nstart_sup
1555 jhpb(nhpb)=j+nstart_seq-nstart_sup
1557 dhpb(nhpb)=dist(i,j)
1560 cd write (iout,'(a)') 'Distance constraints:'
1565 cd if (ii.gt.nres) then
1570 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1571 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1572 cd & dhpb(i),forcon(i)
1576 c----------------------------------------------------------------------------
1578 implicit real*8 (a-h,o-z)
1579 include 'DIMENSIONS'
1580 include 'COMMON.MAP'
1581 include 'COMMON.IOUNITS'
1582 character*3 angid(4) /'THE','PHI','ALP','OME'/
1583 character*80 mapcard,ucase
1585 read (inp,'(a)') mapcard
1586 mapcard=ucase(mapcard)
1587 if (index(mapcard,'PHI').gt.0) then
1589 else if (index(mapcard,'THE').gt.0) then
1591 else if (index(mapcard,'ALP').gt.0) then
1593 else if (index(mapcard,'OME').gt.0) then
1596 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1597 stop 'Error - illegal variable spec in MAP card.'
1599 call readi (mapcard,'RES1',res1(imap),0)
1600 call readi (mapcard,'RES2',res2(imap),0)
1601 if (res1(imap).eq.0) then
1602 res1(imap)=res2(imap)
1603 else if (res2(imap).eq.0) then
1604 res2(imap)=res1(imap)
1606 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1608 & 'Error - illegal definition of variable group in MAP.'
1609 stop 'Error - illegal definition of variable group in MAP.'
1611 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1612 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1613 call readi(mapcard,'NSTEP',nstep(imap),0)
1614 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1616 & 'Illegal boundary and/or step size specification in MAP.'
1617 stop 'Illegal boundary and/or step size specification in MAP.'
1622 c----------------------------------------------------------------------------
1624 implicit real*8 (a-h,o-z)
1625 include 'DIMENSIONS'
1626 include 'COMMON.IOUNITS'
1627 include 'COMMON.GEO'
1628 include 'COMMON.CSA'
1629 include 'COMMON.BANK'
1630 include 'COMMON.CONTROL'
1632 character*620 mcmcard
1633 call card_concat(mcmcard)
1635 call readi(mcmcard,'NCONF',nconf,50)
1636 call readi(mcmcard,'NADD',nadd,0)
1637 call readi(mcmcard,'JSTART',jstart,1)
1638 call readi(mcmcard,'JEND',jend,1)
1639 call readi(mcmcard,'NSTMAX',nstmax,500000)
1640 call readi(mcmcard,'N0',n0,1)
1641 call readi(mcmcard,'N1',n1,6)
1642 call readi(mcmcard,'N2',n2,4)
1643 call readi(mcmcard,'N3',n3,0)
1644 call readi(mcmcard,'N4',n4,0)
1645 call readi(mcmcard,'N5',n5,0)
1646 call readi(mcmcard,'N6',n6,10)
1647 call readi(mcmcard,'N7',n7,0)
1648 call readi(mcmcard,'N8',n8,0)
1649 call readi(mcmcard,'N9',n9,0)
1650 call readi(mcmcard,'N14',n14,0)
1651 call readi(mcmcard,'N15',n15,0)
1652 call readi(mcmcard,'N16',n16,0)
1653 call readi(mcmcard,'N17',n17,0)
1654 call readi(mcmcard,'N18',n18,0)
1656 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1658 call readi(mcmcard,'NDIFF',ndiff,2)
1659 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1660 call readi(mcmcard,'IS1',is1,1)
1661 call readi(mcmcard,'IS2',is2,8)
1662 call readi(mcmcard,'NRAN0',nran0,4)
1663 call readi(mcmcard,'NRAN1',nran1,2)
1664 call readi(mcmcard,'IRR',irr,1)
1665 call readi(mcmcard,'NSEED',nseed,20)
1666 call readi(mcmcard,'NTOTAL',ntotal,10000)
1667 call reada(mcmcard,'CUT1',cut1,2.0d0)
1668 call reada(mcmcard,'CUT2',cut2,5.0d0)
1669 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1670 call readi(mcmcard,'ICMAX',icmax,3)
1671 call readi(mcmcard,'IRESTART',irestart,0)
1672 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1675 call reada(mcmcard,'DELE',dele,20.0d0)
1676 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1677 call readi(mcmcard,'IREF',iref,0)
1678 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1679 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1680 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1681 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1682 write (iout,*) "NCONF_IN",nconf_in
1685 c----------------------------------------------------------------------------
1686 cfmc subroutine mcmfread
1687 cfmc implicit real*8 (a-h,o-z)
1688 cfmc include 'DIMENSIONS'
1689 cfmc include 'COMMON.MCMF'
1690 cfmc include 'COMMON.IOUNITS'
1691 cfmc include 'COMMON.GEO'
1692 cfmc character*80 ucase
1693 cfmc character*620 mcmcard
1694 cfmc call card_concat(mcmcard)
1696 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1697 cfmc write(iout,*)'MAXRANT=',maxrant
1698 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1699 cfmc write(iout,*)'MAXFAM=',maxfam
1700 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1701 cfmc write(iout,*)'NNET1=',nnet1
1702 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1703 cfmc write(iout,*)'NNET2=',nnet2
1704 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1705 cfmc write(iout,*)'NNET3=',nnet3
1706 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1707 cfmc write(iout,*)'ILASTT=',ilastt
1708 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1709 cfmc write(iout,*)'MAXSTR=',maxstr
1710 cfmc maxstr_f=maxstr/maxfam
1711 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1712 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1713 cfmc write(iout,*)'NMCMF=',nmcmf
1714 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1715 cfmc write(iout,*)'IFOCUS=',ifocus
1716 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1717 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1718 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1719 cfmc write(iout,*)'INTPRT=',intprt
1720 cfmc call readi(mcmcard,'IPRT',iprt,100)
1721 cfmc write(iout,*)'IPRT=',iprt
1722 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1723 cfmc write(iout,*)'IMAXTR=',imaxtr
1724 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1725 cfmc write(iout,*)'MAXEVEN=',maxeven
1726 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1727 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1728 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1729 cfmc write(iout,*)'INIMIN=',inimin
1730 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1731 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1732 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1733 cfmc write(iout,*)'NTHREAD=',nthread
1734 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1735 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1736 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1737 cfmc write(iout,*)'MAXPERT=',maxpert
1738 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1739 cfmc write(iout,*)'IRMSD=',irmsd
1740 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1741 cfmc write(iout,*)'DENEMIN=',denemin
1742 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1743 cfmc write(iout,*)'RCUT1S=',rcut1s
1744 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1745 cfmc write(iout,*)'RCUT1E=',rcut1e
1746 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1747 cfmc write(iout,*)'RCUT2S=',rcut2s
1748 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1749 cfmc write(iout,*)'RCUT2E=',rcut2e
1750 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1751 cfmc write(iout,*)'DPERT1=',d_pert1
1752 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1753 cfmc write(iout,*)'DPERT1A=',d_pert1a
1754 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1755 cfmc write(iout,*)'DPERT2=',d_pert2
1756 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1757 cfmc write(iout,*)'DPERT2A=',d_pert2a
1758 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1759 cfmc write(iout,*)'DPERT2B=',d_pert2b
1760 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1761 cfmc write(iout,*)'DPERT2C=',d_pert2c
1762 cfmc d_pert1=deg2rad*d_pert1
1763 cfmc d_pert1a=deg2rad*d_pert1a
1764 cfmc d_pert2=deg2rad*d_pert2
1765 cfmc d_pert2a=deg2rad*d_pert2a
1766 cfmc d_pert2b=deg2rad*d_pert2b
1767 cfmc d_pert2c=deg2rad*d_pert2c
1768 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1769 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1770 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1771 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1772 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1773 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1774 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1775 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1776 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1777 cfmc write(iout,*)'RCUTINI=',rcutini
1778 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1779 cfmc write(iout,*)'GRAT=',grat
1780 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1781 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1785 c----------------------------------------------------------------------------
1787 implicit real*8 (a-h,o-z)
1788 include 'DIMENSIONS'
1789 include 'COMMON.MCM'
1790 include 'COMMON.MCE'
1791 include 'COMMON.IOUNITS'
1793 character*320 mcmcard
1794 call card_concat(mcmcard)
1795 call readi(mcmcard,'MAXACC',maxacc,100)
1796 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1797 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1798 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1799 call readi(mcmcard,'MAXREPM',maxrepm,200)
1800 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1801 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1802 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1803 call reada(mcmcard,'E_UP',e_up,5.0D0)
1804 call reada(mcmcard,'DELTE',delte,0.1D0)
1805 call readi(mcmcard,'NSWEEP',nsweep,5)
1806 call readi(mcmcard,'NSTEPH',nsteph,0)
1807 call readi(mcmcard,'NSTEPC',nstepc,0)
1808 call reada(mcmcard,'TMIN',tmin,298.0D0)
1809 call reada(mcmcard,'TMAX',tmax,298.0D0)
1810 call readi(mcmcard,'NWINDOW',nwindow,0)
1811 call readi(mcmcard,'PRINT_MC',print_mc,0)
1812 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1813 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1814 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1815 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1816 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1817 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1818 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1819 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1820 if (nwindow.gt.0) then
1821 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1823 winlen(i)=winend(i)-winstart(i)+1
1826 if (tmax.lt.tmin) tmax=tmin
1827 if (tmax.eq.tmin) then
1831 if (nstepc.gt.0 .and. nsteph.gt.0) then
1832 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1833 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1835 C Probabilities of different move types
1836 sumpro_type(0)=0.0D0
1837 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1838 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1839 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1840 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1841 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1842 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1843 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1845 print *,'i',i,' sumprotype',sumpro_type(i)
1846 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1847 print *,'i',i,' sumprotype',sumpro_type(i)
1851 c----------------------------------------------------------------------------
1852 subroutine read_minim
1853 implicit real*8 (a-h,o-z)
1854 include 'DIMENSIONS'
1855 include 'COMMON.MINIM'
1856 include 'COMMON.IOUNITS'
1857 include 'COMMON.CONTROL'
1858 include 'COMMON.SETUP'
1860 character*320 minimcard
1861 call card_concat(minimcard)
1862 call readi(minimcard,'MAXMIN',maxmin,2000)
1863 call readi(minimcard,'MAXFUN',maxfun,5000)
1864 call readi(minimcard,'MINMIN',minmin,maxmin)
1865 call readi(minimcard,'MINFUN',minfun,maxmin)
1866 call reada(minimcard,'TOLF',tolf,1.0D-2)
1867 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1868 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1869 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1870 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1872 if (.not. out1file .or. me.eq.king) then
1874 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1875 & 'Options in energy minimization:'
1876 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1877 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1878 & 'MinMin:',MinMin,' MinFun:',MinFun,
1879 & ' TolF:',TolF,' RTolF:',RTolF
1885 c----------------------------------------------------------------------------
1886 subroutine read_angles(kanal,*)
1887 implicit real*8 (a-h,o-z)
1888 include 'DIMENSIONS'
1889 include 'COMMON.GEO'
1890 include 'COMMON.VAR'
1891 include 'COMMON.CHAIN'
1892 include 'COMMON.IOUNITS'
1893 include 'COMMON.CONTROL'
1894 c Read angles from input
1896 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1897 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1898 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1899 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1902 c 9/7/01 avoid 180 deg valence angle
1903 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1905 theta(i)=deg2rad*theta(i)
1906 phi(i)=deg2rad*phi(i)
1907 alph(i)=deg2rad*alph(i)
1908 omeg(i)=deg2rad*omeg(i)
1913 c----------------------------------------------------------------------------
1914 subroutine reada(rekord,lancuch,wartosc,default)
1916 character*(*) rekord,lancuch
1917 double precision wartosc,default
1920 iread=index(rekord,lancuch)
1921 if (iread.eq.0) then
1925 iread=iread+ilen(lancuch)+1
1926 read (rekord(iread:),*,err=10,end=10) wartosc
1931 c----------------------------------------------------------------------------
1932 subroutine readi(rekord,lancuch,wartosc,default)
1934 character*(*) rekord,lancuch
1935 integer wartosc,default
1938 iread=index(rekord,lancuch)
1939 if (iread.eq.0) then
1943 iread=iread+ilen(lancuch)+1
1944 read (rekord(iread:),*,err=10,end=10) wartosc
1949 c----------------------------------------------------------------------------
1950 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1953 integer tablica(dim),default
1954 character*(*) rekord,lancuch
1961 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1962 if (iread.eq.0) return
1963 iread=iread+ilen(lancuch)+1
1964 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1967 c----------------------------------------------------------------------------
1968 subroutine multreada(rekord,lancuch,tablica,dim,default)
1971 double precision tablica(dim),default
1972 character*(*) rekord,lancuch
1979 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1980 if (iread.eq.0) return
1981 iread=iread+ilen(lancuch)+1
1982 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1985 c----------------------------------------------------------------------------
1986 subroutine openunits
1987 implicit real*8 (a-h,o-z)
1988 include 'DIMENSIONS'
1991 character*16 form,nodename
1994 include 'COMMON.SETUP'
1995 include 'COMMON.IOUNITS'
1997 include 'COMMON.CONTROL'
1998 integer lenpre,lenpot,ilen,lentmp
2000 character*3 out1file_text,ucase
2003 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2004 call getenv_loc("PREFIX",prefix)
2006 call getenv_loc("POT",pot)
2007 call getenv_loc("DIRTMP",tmpdir)
2008 call getenv_loc("CURDIR",curdir)
2009 call getenv_loc("OUT1FILE",out1file_text)
2010 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2011 out1file_text=ucase(out1file_text)
2012 if (out1file_text(1:1).eq."Y") then
2015 out1file=fg_rank.gt.0
2020 if (lentmp.gt.0) then
2021 write (*,'(80(1h!))')
2022 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2023 write (*,'(80(1h!))')
2024 write (*,*)"All output files will be on node /tmp directory."
2026 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2027 if (me.eq.king) then
2028 write (*,*) "The master node is ",nodename
2029 else if (fg_rank.eq.0) then
2030 write (*,*) "I am the CG slave node ",nodename
2032 write (*,*) "I am the FG slave node ",nodename
2035 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2036 lenpre = lentmp+lenpre+1
2038 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2039 C Get the names and open the input files
2040 #if defined(WINIFL) || defined(WINPGI)
2041 open(1,file=pref_orig(:ilen(pref_orig))//
2042 & '.inp',status='old',readonly,shared)
2043 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2044 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2045 C Get parameter filenames and open the parameter files.
2046 call getenv_loc('BONDPAR',bondname)
2047 open (ibond,file=bondname,status='old',readonly,shared)
2048 call getenv_loc('THETPAR',thetname)
2049 open (ithep,file=thetname,status='old',readonly,shared)
2051 call getenv_loc('THETPARPDB',thetname_pdb)
2052 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2054 call getenv_loc('ROTPAR',rotname)
2055 open (irotam,file=rotname,status='old',readonly,shared)
2057 call getenv_loc('ROTPARPDB',rotname_pdb)
2058 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2060 call getenv_loc('TORPAR',torname)
2061 open (itorp,file=torname,status='old',readonly,shared)
2062 call getenv_loc('TORDPAR',tordname)
2063 open (itordp,file=tordname,status='old',readonly,shared)
2064 call getenv_loc('FOURIER',fouriername)
2065 open (ifourier,file=fouriername,status='old',readonly,shared)
2066 call getenv_loc('ELEPAR',elename)
2067 open (ielep,file=elename,status='old',readonly,shared)
2068 call getenv_loc('SIDEPAR',sidename)
2069 open (isidep,file=sidename,status='old',readonly,shared)
2070 #elif (defined CRAY) || (defined AIX)
2071 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2073 c print *,"Processor",myrank," opened file 1"
2074 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2075 c print *,"Processor",myrank," opened file 9"
2076 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2077 C Get parameter filenames and open the parameter files.
2078 call getenv_loc('BONDPAR',bondname)
2079 open (ibond,file=bondname,status='old',action='read')
2080 c print *,"Processor",myrank," opened file IBOND"
2081 call getenv_loc('THETPAR',thetname)
2082 open (ithep,file=thetname,status='old',action='read')
2083 c print *,"Processor",myrank," opened file ITHEP"
2085 call getenv_loc('THETPARPDB',thetname_pdb)
2086 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2088 call getenv_loc('ROTPAR',rotname)
2089 open (irotam,file=rotname,status='old',action='read')
2090 c print *,"Processor",myrank," opened file IROTAM"
2092 call getenv_loc('ROTPARPDB',rotname_pdb)
2093 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2095 call getenv_loc('TORPAR',torname)
2096 open (itorp,file=torname,status='old',action='read')
2097 c print *,"Processor",myrank," opened file ITORP"
2098 call getenv_loc('TORDPAR',tordname)
2099 open (itordp,file=tordname,status='old',action='read')
2100 c print *,"Processor",myrank," opened file ITORDP"
2101 call getenv_loc('SCCORPAR',sccorname)
2102 open (isccor,file=sccorname,status='old',action='read')
2103 c print *,"Processor",myrank," opened file ISCCOR"
2104 call getenv_loc('FOURIER',fouriername)
2105 open (ifourier,file=fouriername,status='old',action='read')
2106 c print *,"Processor",myrank," opened file IFOURIER"
2107 call getenv_loc('ELEPAR',elename)
2108 open (ielep,file=elename,status='old',action='read')
2109 c print *,"Processor",myrank," opened file IELEP"
2110 call getenv_loc('SIDEPAR',sidename)
2111 open (isidep,file=sidename,status='old',action='read')
2112 c print *,"Processor",myrank," opened file ISIDEP"
2113 c print *,"Processor",myrank," opened parameter files"
2115 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2116 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2117 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2118 C Get parameter filenames and open the parameter files.
2119 call getenv_loc('BONDPAR',bondname)
2120 open (ibond,file=bondname,status='old')
2121 call getenv_loc('THETPAR',thetname)
2122 open (ithep,file=thetname,status='old')
2124 call getenv_loc('THETPARPDB',thetname_pdb)
2125 open (ithep_pdb,file=thetname_pdb,status='old')
2127 call getenv_loc('ROTPAR',rotname)
2128 open (irotam,file=rotname,status='old')
2130 call getenv_loc('ROTPARPDB',rotname_pdb)
2131 open (irotam_pdb,file=rotname_pdb,status='old')
2133 call getenv_loc('TORPAR',torname)
2134 open (itorp,file=torname,status='old')
2135 call getenv_loc('TORDPAR',tordname)
2136 open (itordp,file=tordname,status='old')
2137 call getenv_loc('SCCORPAR',sccorname)
2138 open (isccor,file=sccorname,status='old')
2139 call getenv_loc('FOURIER',fouriername)
2140 open (ifourier,file=fouriername,status='old')
2141 call getenv_loc('ELEPAR',elename)
2142 open (ielep,file=elename,status='old')
2143 call getenv_loc('SIDEPAR',sidename)
2144 open (isidep,file=sidename,status='old')
2145 call getenv_loc('LIPTRANPAR',liptranname)
2146 open (iliptranpar,file=liptranname,status='old')
2148 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2150 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2151 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2152 C Get parameter filenames and open the parameter files.
2153 call getenv_loc('BONDPAR',bondname)
2154 open (ibond,file=bondname,status='old',readonly)
2155 call getenv_loc('THETPAR',thetname)
2156 open (ithep,file=thetname,status='old',readonly)
2157 call getenv_loc('ROTPAR',rotname)
2158 open (irotam,file=rotname,status='old',readonly)
2159 call getenv_loc('TORPAR',torname)
2160 open (itorp,file=torname,status='old',readonly)
2161 call getenv_loc('TORDPAR',tordname)
2162 open (itordp,file=tordname,status='old',readonly)
2163 call getenv_loc('SCCORPAR',sccorname)
2164 open (isccor,file=sccorname,status='old',readonly)
2166 call getenv_loc('THETPARPDB',thetname_pdb)
2167 c print *,"thetname_pdb ",thetname_pdb
2168 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2169 c print *,ithep_pdb," opened"
2171 call getenv_loc('FOURIER',fouriername)
2172 open (ifourier,file=fouriername,status='old',readonly)
2173 call getenv_loc('ELEPAR',elename)
2174 open (ielep,file=elename,status='old',readonly)
2175 call getenv_loc('SIDEPAR',sidename)
2176 open (isidep,file=sidename,status='old',readonly)
2177 call getenv_loc('LIPTRANPAR',liptranname)
2178 open (iliptranpar,file=liptranname,status='old',readonly)
2180 call getenv_loc('ROTPARPDB',rotname_pdb)
2181 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2186 C 8/9/01 In the newest version SCp interaction constants are read from a file
2187 C Use -DOLDSCP to use hard-coded constants instead.
2189 call getenv_loc('SCPPAR',scpname)
2190 #if defined(WINIFL) || defined(WINPGI)
2191 open (iscpp,file=scpname,status='old',readonly,shared)
2192 #elif (defined CRAY) || (defined AIX)
2193 open (iscpp,file=scpname,status='old',action='read')
2195 open (iscpp,file=scpname,status='old')
2197 open (iscpp,file=scpname,status='old',readonly)
2200 call getenv_loc('PATTERN',patname)
2201 #if defined(WINIFL) || defined(WINPGI)
2202 open (icbase,file=patname,status='old',readonly,shared)
2203 #elif (defined CRAY) || (defined AIX)
2204 open (icbase,file=patname,status='old',action='read')
2206 open (icbase,file=patname,status='old')
2208 open (icbase,file=patname,status='old',readonly)
2211 C Open output file only for CG processes
2212 c print *,"Processor",myrank," fg_rank",fg_rank
2213 if (fg_rank.eq.0) then
2215 if (nodes.eq.1) then
2218 npos = dlog10(dfloat(nodes-1))+1
2220 if (npos.lt.3) npos=3
2221 write (liczba,'(i1)') npos
2222 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2224 write (liczba,form) me
2225 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2226 & liczba(:ilen(liczba))
2227 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2229 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2231 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2232 & liczba(:ilen(liczba))//'.mol2'
2233 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2234 & liczba(:ilen(liczba))//'.stat'
2236 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2237 & //liczba(:ilen(liczba))//'.stat')
2238 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2241 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2242 & liczba(:ilen(liczba))//'.const'
2247 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2248 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2249 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2250 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2251 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2253 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2255 rest2name=prefix(:ilen(prefix))//'.rst'
2257 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2260 #if defined(AIX) || defined(PGI)
2261 if (me.eq.king .or. .not. out1file)
2262 & open(iout,file=outname,status='unknown')
2264 if (fg_rank.gt.0) then
2265 write (liczba,'(i3.3)') myrank/nfgtasks
2266 write (ll,'(bz,i3.3)') fg_rank
2267 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2272 open(igeom,file=intname,status='unknown',position='append')
2273 open(ipdb,file=pdbname,status='unknown')
2274 open(imol2,file=mol2name,status='unknown')
2275 open(istat,file=statname,status='unknown',position='append')
2277 c1out open(iout,file=outname,status='unknown')
2280 if (me.eq.king .or. .not.out1file)
2281 & open(iout,file=outname,status='unknown')
2283 if (fg_rank.gt.0) then
2284 write (liczba,'(i3.3)') myrank/nfgtasks
2285 write (ll,'(bz,i3.3)') fg_rank
2286 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2291 open(igeom,file=intname,status='unknown',access='append')
2292 open(ipdb,file=pdbname,status='unknown')
2293 open(imol2,file=mol2name,status='unknown')
2294 open(istat,file=statname,status='unknown',access='append')
2296 c1out open(iout,file=outname,status='unknown')
2299 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2300 csa_seed=prefix(:lenpre)//'.CSA.seed'
2301 csa_history=prefix(:lenpre)//'.CSA.history'
2302 csa_bank=prefix(:lenpre)//'.CSA.bank'
2303 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2304 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2305 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2306 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2307 csa_int=prefix(:lenpre)//'.int'
2308 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2309 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2310 csa_in=prefix(:lenpre)//'.CSA.in'
2311 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2314 write (iout,'(80(1h-))')
2315 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2316 write (iout,'(80(1h-))')
2317 write (iout,*) "Input file : ",
2318 & pref_orig(:ilen(pref_orig))//'.inp'
2319 write (iout,*) "Output file : ",
2320 & outname(:ilen(outname))
2322 write (iout,*) "Sidechain potential file : ",
2323 & sidename(:ilen(sidename))
2325 write (iout,*) "SCp potential file : ",
2326 & scpname(:ilen(scpname))
2328 write (iout,*) "Electrostatic potential file : ",
2329 & elename(:ilen(elename))
2330 write (iout,*) "Cumulant coefficient file : ",
2331 & fouriername(:ilen(fouriername))
2332 write (iout,*) "Torsional parameter file : ",
2333 & torname(:ilen(torname))
2334 write (iout,*) "Double torsional parameter file : ",
2335 & tordname(:ilen(tordname))
2336 write (iout,*) "SCCOR parameter file : ",
2337 & sccorname(:ilen(sccorname))
2338 write (iout,*) "Bond & inertia constant file : ",
2339 & bondname(:ilen(bondname))
2340 write (iout,*) "Bending parameter file : ",
2341 & thetname(:ilen(thetname))
2342 write (iout,*) "Rotamer parameter file : ",
2343 & rotname(:ilen(rotname))
2344 write (iout,*) "Thetpdb parameter file : ",
2345 & thetname_pdb(:ilen(thetname_pdb))
2346 write (iout,*) "Threading database : ",
2347 & patname(:ilen(patname))
2349 &write (iout,*)" DIRTMP : ",
2351 write (iout,'(80(1h-))')
2355 c----------------------------------------------------------------------------
2356 subroutine card_concat(card)
2357 implicit real*8 (a-h,o-z)
2358 include 'DIMENSIONS'
2359 include 'COMMON.IOUNITS'
2361 character*80 karta,ucase
2363 read (inp,'(a)') karta
2366 do while (karta(80:80).eq.'&')
2367 card=card(:ilen(card)+1)//karta(:79)
2368 read (inp,'(a)') karta
2371 card=card(:ilen(card)+1)//karta
2374 c----------------------------------------------------------------------------------
2376 implicit real*8 (a-h,o-z)
2377 include 'DIMENSIONS'
2378 include 'COMMON.CHAIN'
2379 include 'COMMON.IOUNITS'
2381 include 'COMMON.CONTROL'
2382 open(irest2,file=rest2name,status='unknown')
2383 read(irest2,*) totT,EK,potE,totE,t_bath
2386 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2389 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2391 if(usampl.or.homol_nset.gt.1) then
2392 read (irest2,*) iset
2397 c---------------------------------------------------------------------------------
2398 subroutine read_fragments
2399 implicit real*8 (a-h,o-z)
2400 include 'DIMENSIONS'
2404 include 'COMMON.SETUP'
2405 include 'COMMON.CHAIN'
2406 include 'COMMON.IOUNITS'
2408 include 'COMMON.CONTROL'
2410 read(inp,*) nset,nfrag,npair,nfrag_back
2411 if(me.eq.king.or..not.out1file)
2412 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2413 & " nfrag_back",nfrag_back
2415 read(inp,*) mset(iset1)
2417 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2419 if(me.eq.king.or..not.out1file)
2420 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2421 & ifrag(2,i,iset1), qinfrag(i,iset1)
2424 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2426 if(me.eq.king.or..not.out1file)
2427 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2428 & ipair(2,i,iset1), qinpair(i,iset1)
2431 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2432 & wfrag_back(3,i,iset1),
2433 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2434 if(me.eq.king.or..not.out1file)
2435 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2436 & wfrag_back(2,i,iset1),
2437 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2438 & ifrag_back(2,i,iset1)
2443 C---------------------------------------------------------------------------
2444 subroutine read_afminp
2445 implicit real*8 (a-h,o-z)
2446 include 'DIMENSIONS'
2450 include 'COMMON.SETUP'
2451 include 'COMMON.CONTROL'
2452 include 'COMMON.CHAIN'
2453 include 'COMMON.IOUNITS'
2454 include 'COMMON.SBRIDGE'
2455 character*320 afmcard
2456 c print *, "wchodze"
2457 call card_concat(afmcard)
2458 call readi(afmcard,"BEG",afmbeg,0)
2459 call readi(afmcard,"END",afmend,0)
2460 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2461 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2462 print *,'FORCE=' ,forceAFMconst
2463 CCCC NOW PROPERTIES FOR AFM
2466 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2468 distafminit=dsqrt(distafminit)
2469 print *,'initdist',distafminit
2473 c-------------------------------------------------------------------------------
2474 subroutine read_saxs_constr
2475 implicit real*8 (a-h,o-z)
2476 include 'DIMENSIONS'
2480 include 'COMMON.SETUP'
2481 include 'COMMON.CONTROL'
2482 include 'COMMON.CHAIN'
2483 include 'COMMON.IOUNITS'
2484 include 'COMMON.SBRIDGE'
2485 double precision cm(3)
2487 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2489 if (saxs_mode.eq.0) then
2490 c SAXS distance distribution
2492 read(inp,*) distsaxs(i),Psaxs(i)
2496 Cnorm = Cnorm + Psaxs(i)
2498 write (iout,*) "Cnorm",Cnorm
2500 Psaxs(i)=Psaxs(i)/Cnorm
2502 write (iout,*) "Normalized distance distribution from SAXS"
2504 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2508 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2510 write (iout,*) "Wsaxs0",Wsaxs0
2514 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2521 cm(j)=cm(j)+Csaxs(j,i)
2529 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2532 write (iout,*) "SAXS sphere coordinates"
2534 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2539 c-------------------------------------------------------------------------------
2540 subroutine read_dist_constr
2541 implicit real*8 (a-h,o-z)
2542 include 'DIMENSIONS'
2546 include 'COMMON.SETUP'
2547 include 'COMMON.CONTROL'
2548 include 'COMMON.CHAIN'
2549 include 'COMMON.IOUNITS'
2550 include 'COMMON.SBRIDGE'
2551 integer ifrag_(2,100),ipair_(2,100)
2552 double precision wfrag_(100),wpair_(100)
2553 character*500 controlcard
2554 logical normalize,next
2556 double precision xlink(4,0:4) /
2558 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2559 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2560 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2561 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2562 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2563 c print *, "WCHODZE"
2564 write (iout,*) "Calling read_dist_constr"
2565 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2571 call card_concat(controlcard)
2572 next = index(controlcard,"NEXT").gt.0
2573 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2574 write (iout,*) "restr_type",restr_type
2575 call readi(controlcard,"NFRAG",nfrag_,0)
2576 call readi(controlcard,"NFRAG",nfrag_,0)
2577 call readi(controlcard,"NPAIR",npair_,0)
2578 call readi(controlcard,"NDIST",ndist_,0)
2579 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2580 if (restr_type.eq.10)
2581 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2582 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2583 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2584 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2585 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2586 normalize = index(controlcard,"NORMALIZE").gt.0
2587 write (iout,*) "WBOLTZD",wboltzd
2588 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2589 c write (iout,*) "IFRAG"
2591 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2593 c write (iout,*) "IPAIR"
2595 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2597 if (nfrag_.gt.0) write (iout,*)
2598 & "Distance restraints as generated from reference structure"
2600 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2601 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2602 & ifrag_(2,i)=nstart_sup+nsup-1
2603 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2605 if (wfrag_(i).eq.0.0d0) cycle
2606 do j=ifrag_(1,i),ifrag_(2,i)-1
2607 do k=j+1,ifrag_(2,i)
2608 c write (iout,*) "j",j," k",k
2610 if (restr_type.eq.1) then
2616 forcon(nhpb)=wfrag_(i)
2617 else if (constr_dist.eq.2) then
2618 if (ddjk.le.dist_cut) then
2624 forcon(nhpb)=wfrag_(i)
2626 else if (restr_type.eq.3) then
2632 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2635 if (.not.out1file .or. me.eq.king)
2636 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2637 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2639 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2640 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2646 if (wpair_(i).eq.0.0d0) cycle
2654 do j=ifrag_(1,ii),ifrag_(2,ii)
2655 do k=ifrag_(1,jj),ifrag_(2,jj)
2656 if (restr_type.eq.1) then
2662 forcon(nhpb)=wfrag_(i)
2663 else if (constr_dist.eq.2) then
2664 if (ddjk.le.dist_cut) then
2670 forcon(nhpb)=wfrag_(i)
2672 else if (restr_type.eq.3) then
2678 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2681 if (.not.out1file .or. me.eq.king)
2682 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
2683 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2685 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
2686 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2693 write (iout,*) "Distance restraints as read from input"
2695 if (restr_type.eq.11) then
2696 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2697 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2698 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2699 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2701 irestr_type(nhpb)=11
2703 if (.not.out1file .or. me.eq.king)
2704 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2705 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2706 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2708 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2709 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2710 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2712 if (ibecarb(nhpb).gt.0) then
2713 ihpb(nhpb)=ihpb(nhpb)+nres
2714 jhpb(nhpb)=jhpb(nhpb)+nres
2716 else if (constr_dist.eq.10) then
2717 c Cross-lonk Markov-like potential
2718 call card_concat(controlcard)
2719 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2720 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2722 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2723 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2724 if (index(controlcard,"ZL").gt.0) then
2726 else if (index(controlcard,"ADH").gt.0) then
2728 else if (index(controlcard,"PDH").gt.0) then
2730 else if (index(controlcard,"DSS").gt.0) then
2735 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2736 & xlink(1,link_type))
2737 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2738 & xlink(2,link_type))
2739 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2740 & xlink(3,link_type))
2741 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2742 & xlink(4,link_type))
2743 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2744 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2745 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2746 if (forcon(nhpb+1).le.0.0d0 .or.
2747 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2749 irestr_type(nhpb)=10
2750 if (ibecarb(nhpb).gt.0) then
2751 ihpb(nhpb)=ihpb(nhpb)+nres
2752 jhpb(nhpb)=jhpb(nhpb)+nres
2755 if (.not.out1file .or. me.eq.king)
2756 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2757 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2758 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2761 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2762 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2763 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2768 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2769 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2770 if (forcon(nhpb+1).gt.0.0d0) then
2772 if (dhpb1(nhpb).eq.0.0d0) then
2777 if (ibecarb(nhpb).gt.0) then
2778 ihpb(nhpb)=ihpb(nhpb)+nres
2779 jhpb(nhpb)=jhpb(nhpb)+nres
2781 if (dhpb(nhpb).eq.0.0d0)
2782 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2785 if (.not.out1file .or. me.eq.king)
2786 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2787 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2789 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2790 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2793 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2794 C if (forcon(nhpb+1).gt.0.0d0) then
2796 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2804 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2805 & fordepthmax=fordepth(i)
2808 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2813 c-------------------------------------------------------------------------------
2815 subroutine read_constr_homology
2817 include 'DIMENSIONS'
2821 include 'COMMON.SETUP'
2822 include 'COMMON.CONTROL'
2823 include 'COMMON.CHAIN'
2824 include 'COMMON.IOUNITS'
2826 include 'COMMON.GEO'
2827 include 'COMMON.INTERACT'
2828 include 'COMMON.NAMES'
2830 c For new homol impl
2832 include 'COMMON.VAR'
2835 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2837 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2838 c & sigma_odl_temp(maxres,maxres,max_template)
2840 character*24 model_ki_dist, model_ki_angle
2841 character*500 controlcard
2842 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2843 logical lprn /.true./
2848 c FP - Nov. 2014 Temporary specifications for new vars
2850 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2852 double precision, dimension (max_template,maxres) :: rescore
2853 double precision, dimension (max_template,maxres) :: rescore2
2854 double precision, dimension (max_template,maxres) :: rescore3
2855 character*24 pdbfile,tpl_k_rescore
2856 c -----------------------------------------------------------------
2857 c Reading multiple PDB ref structures and calculation of retraints
2858 c not using pre-computed ones stored in files model_ki_{dist,angle}
2860 c -----------------------------------------------------------------
2863 c Alternative: reading from input
2864 call card_concat(controlcard)
2865 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2866 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2867 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2868 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2869 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2870 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2871 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2872 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2873 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2874 if(.not.read2sigma.and.start_from_model) then
2875 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2876 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2877 start_from_model=.false.
2879 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2880 & write(iout,*) 'START_FROM_MODELS is ON'
2881 if(start_from_model .and. rest) then
2882 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2883 write(iout,*) 'START_FROM_MODELS is OFF'
2884 write(iout,*) 'remove restart keyword from input'
2887 if (homol_nset.gt.1)then
2888 call card_concat(controlcard)
2889 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2890 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2891 write(iout,*) "iset homology_weight "
2893 write(iout,*) i,waga_homology(i)
2896 iset=mod(kolor,homol_nset)+1
2899 waga_homology(1)=1.0
2902 cd write (iout,*) "nnt",nnt," nct",nct
2909 c write(iout,*) 'nnt=',nnt,'nct=',nct
2912 do k=1,constr_homology
2925 if (read_homol_frag) then
2926 call read_klapaucjusz
2929 do k=1,constr_homology
2931 read(inp,'(a)') pdbfile
2932 if(me.eq.king .or. .not. out1file)
2933 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2934 & pdbfile(:ilen(pdbfile))
2935 open(ipdbin,file=pdbfile,status='old',err=33)
2937 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2938 & pdbfile(:ilen(pdbfile))
2941 c print *,'Begin reading pdb data'
2943 c Files containing res sim or local scores (former containing sigmas)
2946 write(kic2,'(bz,i2.2)') k
2948 tpl_k_rescore="template"//kic2//".sco"
2951 if (read2sigma) then
2952 call readpdb_template(k)
2957 c Distance restraints
2960 C Copy the coordinates from reference coordinates (?)
2964 c write (iout,*) "c(",j,i,") =",c(j,i)
2968 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2970 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2971 open (ientin,file=tpl_k_rescore,status='old')
2972 if (nnt.gt.1) rescore(k,1)=0.0d0
2973 do irec=nnt,nct ! loop for reading res sim
2974 if (read2sigma) then
2975 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2976 & rescore3_tmp,idomain_tmp
2978 idomain(k,i_tmp)=idomain_tmp
2979 rescore(k,i_tmp)=rescore_tmp
2980 rescore2(k,i_tmp)=rescore2_tmp
2981 rescore3(k,i_tmp)=rescore3_tmp
2982 if (.not. out1file .or. me.eq.king)
2983 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
2984 & i_tmp,rescore2_tmp,rescore_tmp,
2985 & rescore3_tmp,idomain_tmp
2988 read (ientin,*,end=1401) rescore_tmp
2990 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2991 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2992 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2997 if (waga_dist.ne.0.0d0) then
3005 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3006 c write (iout,*) k,i,j,distal,dist2_cut
3008 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3009 & .and. distal.le.dist2_cut ) then
3015 c write (iout,*) "k",k
3016 c write (iout,*) "i",i," j",j," constr_homology",
3021 if (read2sigma) then
3024 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3026 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3027 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3028 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3030 if (odl(k,ii).le.dist_cut) then
3031 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3034 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3035 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3037 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3038 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3042 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3045 l_homo(k,ii)=.false.
3052 c Theta, dihedral and SC retraints
3054 if (waga_angle.gt.0.0d0) then
3055 c open (ientin,file=tpl_k_sigma_dih,status='old')
3056 c do irec=1,maxres-3 ! loop for reading sigma_dih
3057 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3058 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3059 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3060 c & sigma_dih(k,i+nnt-1)
3065 if (idomain(k,i).eq.0) then
3069 dih(k,i)=phiref(i) ! right?
3070 c read (ientin,*) sigma_dih(k,i) ! original variant
3071 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3072 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3073 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3074 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3075 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3077 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3078 & rescore(k,i-2)+rescore(k,i-3))/4.0
3079 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3080 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3081 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3082 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3083 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3084 c Instead of res sim other local measure of b/b str reliability possible
3085 if (sigma_dih(k,i).ne.0)
3086 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3087 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3092 if (waga_theta.gt.0.0d0) then
3093 c open (ientin,file=tpl_k_sigma_theta,status='old')
3094 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3095 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3096 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3097 c & sigma_theta(k,i+nnt-1)
3102 do i = nnt+2,nct ! right? without parallel.
3103 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3104 c do i=ithet_start,ithet_end ! with FG parallel.
3105 if (idomain(k,i).eq.0) then
3106 sigma_theta(k,i)=0.0
3109 thetatpl(k,i)=thetaref(i)
3110 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3111 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3112 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3113 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3114 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3115 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3116 & rescore(k,i-2))/3.0
3117 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3118 if (sigma_theta(k,i).ne.0)
3119 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3121 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3122 c rescore(k,i-2) ! right expression ?
3123 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3127 if (waga_d.gt.0.0d0) then
3128 c open (ientin,file=tpl_k_sigma_d,status='old')
3129 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3130 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3131 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3132 c & sigma_d(k,i+nnt-1)
3136 do i = nnt,nct ! right? without parallel.
3137 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3138 c do i=loc_start,loc_end ! with FG parallel.
3139 if (itype(i).eq.10) cycle
3140 if (idomain(k,i).eq.0 ) then
3147 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3148 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3149 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3150 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3151 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3152 if (sigma_d(k,i).ne.0)
3153 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3155 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3156 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3157 c read (ientin,*) sigma_d(k,i) ! 1st variant
3162 c remove distance restraints not used in any model from the list
3163 c shift data in all arrays
3165 if (waga_dist.ne.0.0d0) then
3171 if (ii_in_use(ii).eq.0.and.liiflag) then
3175 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3176 & .not.liiflag.and.ii.eq.lim_odl) then
3177 if (ii.eq.lim_odl) then
3178 iishift=ii-iistart+1
3183 do ki=iistart,lim_odl-iishift
3184 ires_homo(ki)=ires_homo(ki+iishift)
3185 jres_homo(ki)=jres_homo(ki+iishift)
3186 ii_in_use(ki)=ii_in_use(ki+iishift)
3187 do k=1,constr_homology
3188 odl(k,ki)=odl(k,ki+iishift)
3189 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3190 l_homo(k,ki)=l_homo(k,ki+iishift)
3194 lim_odl=lim_odl-iishift
3200 endif ! .not. klapaucjusz
3202 if (constr_homology.gt.0) call homology_partition
3203 if (constr_homology.gt.0) call init_int_table
3204 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3205 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3209 if (.not.lprn) return
3210 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3211 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3212 write (iout,*) "Distance restraints from templates"
3214 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3215 & ii,ires_homo(ii),jres_homo(ii),
3216 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3217 & ki=1,constr_homology)
3219 write (iout,*) "Dihedral angle restraints from templates"
3221 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3222 & (rad2deg*dih(ki,i),
3223 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3225 write (iout,*) "Virtual-bond angle restraints from templates"
3227 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3228 & (rad2deg*thetatpl(ki,i),
3229 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3231 write (iout,*) "SC restraints from templates"
3233 write(iout,'(i5,100(4f8.2,4x))') i,
3234 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3235 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3238 c -----------------------------------------------------------------
3241 c----------------------------------------------------------------------
3244 subroutine flush(iu)
3249 subroutine flush(iu)
3254 c------------------------------------------------------------------------------
3255 subroutine copy_to_tmp(source)
3256 include "DIMENSIONS"
3257 include "COMMON.IOUNITS"
3258 character*(*) source
3259 character* 256 tmpfile
3263 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3264 inquire(file=tmpfile,exist=ex)
3266 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3267 & " to temporary directory..."
3268 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3269 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3273 c------------------------------------------------------------------------------
3274 subroutine move_from_tmp(source)
3275 include "DIMENSIONS"
3276 include "COMMON.IOUNITS"
3277 character*(*) source
3280 write (*,*) "Moving ",source(:ilen(source)),
3281 & " from temporary directory to working directory"
3282 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3283 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3286 c------------------------------------------------------------------------------
3287 subroutine random_init(seed)
3289 C Initialize random number generator
3291 implicit real*8 (a-h,o-z)
3292 include 'DIMENSIONS'
3298 logical OKRandom, prng_restart
3300 integer iseed_array(4)
3302 include 'COMMON.IOUNITS'
3303 include 'COMMON.TIME1'
3304 include 'COMMON.THREAD'
3305 include 'COMMON.SBRIDGE'
3306 include 'COMMON.CONTROL'
3307 include 'COMMON.MCM'
3308 include 'COMMON.MAP'
3309 include 'COMMON.HEADER'
3310 include 'COMMON.CSA'
3311 include 'COMMON.CHAIN'
3312 include 'COMMON.MUCA'
3314 include 'COMMON.FFIELD'
3315 include 'COMMON.SETUP'
3316 iseed=-dint(dabs(seed))
3317 if (iseed.eq.0) then
3318 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3319 & 'Random seed undefined. The program will stop.'
3320 write (*,'(/80(1h*)/20x,a/80(1h*))')
3321 & 'Random seed undefined. The program will stop.'
3323 call mpi_finalize(mpi_comm_world,ierr)
3325 stop 'Bad random seed.'
3328 if (fg_rank.eq.0) then
3332 if(me.eq.king .or. .not. out1file)
3333 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3334 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3335 OKRandom = prng_restart(me,iseedi8)
3338 tmp=65536.0d0**(4-i)
3339 iseed_array(i) = dint(seed/tmp)
3340 seed=seed-iseed_array(i)*tmp
3343 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3344 & (iseed_array(i),i=1,4)
3345 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3346 & (iseed_array(i),i=1,4)
3347 OKRandom = prng_restart(me,iseed_array)
3350 c r1 = prng_next(me)
3351 r1=ran_number(0.0D0,1.0D0)
3352 if(me.eq.king .or. .not. out1file)
3353 & write (iout,*) 'ran_num',r1
3354 if (r1.lt.0.0d0) OKRandom=.false.
3356 if (.not.OKRandom) then
3357 write (iout,*) 'PRNG IS NOT WORKING!!!'
3358 print *,'PRNG IS NOT WORKING!!!'
3361 call mpi_abort(mpi_comm_world,error_msg,ierr)
3364 write (iout,*) 'too many processors for parallel prng'
3365 write (*,*) 'too many processors for parallel prng'
3373 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3378 c -----------------------------------------------------------------
3379 subroutine read_klapaucjusz
3381 include 'DIMENSIONS'
3385 include 'COMMON.SETUP'
3386 include 'COMMON.CONTROL'
3387 include 'COMMON.CHAIN'
3388 include 'COMMON.IOUNITS'
3390 include 'COMMON.GEO'
3391 include 'COMMON.INTERACT'
3392 include 'COMMON.NAMES'
3393 character*256 fragfile
3394 integer ninclust(maxclust),inclust(max_template,maxclust),
3395 & nresclust(maxclust),iresclust(maxres,maxclust)
3398 character*24 model_ki_dist, model_ki_angle
3399 character*500 controlcard
3400 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3401 logical lprn /.true./
3407 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3408 double precision, dimension (max_template,maxres) :: rescore
3409 double precision, dimension (max_template,maxres) :: rescore2
3410 character*24 pdbfile,tpl_k_rescore
3413 c For new homol impl
3415 include 'COMMON.VAR'
3417 call getenv("FRAGFILE",fragfile)
3418 open(ientin,file=fragfile,status="old",err=10)
3419 read(ientin,*) constr_homology,nclust
3425 do k=1,constr_homology
3426 read(ientin,'(a)') pdbfile
3427 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3428 & pdbfile(:ilen(pdbfile))
3429 open(ipdbin,file=pdbfile,status='old',err=33)
3431 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3432 & pdbfile(:ilen(pdbfile))
3436 call readpdb_template(k)
3444 read(ientin,*) ninclust(i),nresclust(i)
3445 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3446 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3449 c Loop over clusters
3452 do ll = 1,ninclust(l)
3460 idomain(k,iresclust(i,l)+1) = 1
3462 idomain(k,iresclust(i,l)) = 1
3466 c Distance restraints
3469 C Copy the coordinates from reference coordinates (?)
3473 c write (iout,*) "c(",j,i,") =",c(j,i)
3476 call int_from_cart(.true.,.false.)
3477 call sc_loc_geom(.false.)
3479 thetaref(i)=theta(i)
3482 if (waga_dist.ne.0.0d0) then
3490 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3491 c write (iout,*) k,i,j,distal,dist2_cut
3493 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3494 & .and. distal.le.dist2_cut ) then
3500 c write (iout,*) "k",k
3501 c write (iout,*) "i",i," j",j," constr_homology",
3506 if (read2sigma) then
3509 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3511 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3512 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3513 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3515 if (odl(k,ii).le.dist_cut) then
3516 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3519 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3520 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3522 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3523 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3527 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3530 c l_homo(k,ii)=.false.
3537 c Theta, dihedral and SC retraints
3539 if (waga_angle.gt.0.0d0) then
3541 if (idomain(k,i).eq.0) then
3542 c sigma_dih(k,i)=0.0
3546 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3547 & rescore(k,i-2)+rescore(k,i-3))/4.0
3548 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3549 c & " sigma_dihed",sigma_dih(k,i)
3550 if (sigma_dih(k,i).ne.0)
3551 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3556 if (waga_theta.gt.0.0d0) then
3558 if (idomain(k,i).eq.0) then
3559 c sigma_theta(k,i)=0.0
3562 thetatpl(k,i)=thetaref(i)
3563 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3564 & rescore(k,i-2))/3.0
3565 if (sigma_theta(k,i).ne.0)
3566 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3570 if (waga_d.gt.0.0d0) then
3572 if (itype(i).eq.10) cycle
3573 if (idomain(k,i).eq.0 ) then
3580 sigma_d(k,i)=rescore(k,i)
3581 if (sigma_d(k,i).ne.0)
3582 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3583 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3589 c remove distance restraints not used in any model from the list
3590 c shift data in all arrays
3592 if (waga_dist.ne.0.0d0) then
3598 if (ii_in_use(ii).eq.0.and.liiflag) then
3602 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3603 & .not.liiflag.and.ii.eq.lim_odl) then
3604 if (ii.eq.lim_odl) then
3605 iishift=ii-iistart+1
3610 do ki=iistart,lim_odl-iishift
3611 ires_homo(ki)=ires_homo(ki+iishift)
3612 jres_homo(ki)=jres_homo(ki+iishift)
3613 ii_in_use(ki)=ii_in_use(ki+iishift)
3614 do k=1,constr_homology
3615 odl(k,ki)=odl(k,ki+iishift)
3616 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3617 l_homo(k,ki)=l_homo(k,ki+iishift)
3621 lim_odl=lim_odl-iishift
3628 10 stop "Error infragment file"
3630 c----------------------------------------------------------------------