8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
743 integer :: maxinter,junk,kk,ii,ncatprotparm
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2
750 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
751 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
752 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 ! For printing parameters after they are read set the following in the UNRES
757 ! setenv PRINT_PARM YES
759 ! To print parameters in LaTeX format rather than as ASCII tables:
763 call getenv_loc("PRINT_PARM",lancuch)
764 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
765 call getenv_loc("LATEX",lancuch)
766 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
768 dwa16=2.0d0**(1.0d0/6.0d0)
770 ! Assign virtual-bond length
773 vblinv2=vblinv*vblinv
775 ! Read the virtual-bond parameters, masses, and moments of inertia
776 ! and Stokes' radii of the peptide group and side chains
778 allocate(dsc(ntyp1)) !(ntyp1)
779 allocate(dsc_inv(ntyp1)) !(ntyp1)
780 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
781 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
783 allocate(nbondterm(ntyp)) !(ntyp)
784 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
786 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(long_r_sidechain(ntyp))
788 allocate(short_r_sidechain(ntyp))
793 allocate(msc(ntyp+1)) !(ntyp+1)
794 allocate(isc(ntyp+1)) !(ntyp+1)
795 allocate(restok(ntyp+1)) !(ntyp+1)
797 read (ibond,*) vbldp0,akp,mp,ip,pstok
800 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
801 dsc(i) = vbldsc0(1,i)
805 dsc_inv(i)=1.0D0/dsc(i)
814 allocate(msc(ntyp+1,5)) !(ntyp+1)
815 allocate(isc(ntyp+1,5)) !(ntyp+1)
816 allocate(restok(ntyp+1,5)) !(ntyp+1)
818 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
820 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
821 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
822 dsc(i) = vbldsc0(1,i)
826 dsc_inv(i)=1.0D0/dsc(i)
830 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
833 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
834 ! dsc(i) = vbldsc0_nucl(1,i)
838 ! dsc_inv(i)=1.0D0/dsc(i)
841 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
842 ! do i=1,ntyp_molec(2)
843 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
844 ! aksc_nucl(j,i),abond0_nucl(j,i),&
845 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
846 ! dsc(i) = vbldsc0(1,i)
850 ! dsc_inv(i)=1.0D0/dsc(i)
855 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
856 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
858 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
860 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
861 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
863 write (iout,'(13x,3f10.5)') &
864 vbldsc0(j,i),aksc(j,i),abond0(j,i)
869 read(iion,*) msc(i,5),restok(i,5)
870 print *,msc(i,5),restok(i,5)
874 read (iion,*) ncatprotparm
875 allocate(catprm(ncatprotparm,4))
877 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
880 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
881 !----------------------------------------------------
882 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
883 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
884 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
885 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
886 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
887 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
897 allocate(liptranene(ntyp))
898 !C reading lipid parameters
899 write (iout,*) "iliptranpar",iliptranpar
901 read(iliptranpar,*) pepliptran
904 read(iliptranpar,*) liptranene(i)
905 print *,liptranene(i)
911 ! Read the parameters of the probability distribution/energy expression
912 ! of the virtual-bond valence angles theta
915 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
916 (bthet(j,i,1,1),j=1,2)
917 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
918 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
919 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
923 athet(1,i,1,-1)=athet(1,i,1,1)
924 athet(2,i,1,-1)=athet(2,i,1,1)
925 bthet(1,i,1,-1)=-bthet(1,i,1,1)
926 bthet(2,i,1,-1)=-bthet(2,i,1,1)
927 athet(1,i,-1,1)=-athet(1,i,1,1)
928 athet(2,i,-1,1)=-athet(2,i,1,1)
929 bthet(1,i,-1,1)=bthet(1,i,1,1)
930 bthet(2,i,-1,1)=bthet(2,i,1,1)
934 athet(1,i,-1,-1)=athet(1,-i,1,1)
935 athet(2,i,-1,-1)=-athet(2,-i,1,1)
936 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
937 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
938 athet(1,i,-1,1)=athet(1,-i,1,1)
939 athet(2,i,-1,1)=-athet(2,-i,1,1)
940 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
941 bthet(2,i,-1,1)=bthet(2,-i,1,1)
942 athet(1,i,1,-1)=-athet(1,-i,1,1)
943 athet(2,i,1,-1)=athet(2,-i,1,1)
944 bthet(1,i,1,-1)=bthet(1,-i,1,1)
945 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
950 polthet(j,i)=polthet(j,-i)
953 gthet(j,i)=gthet(j,-i)
961 'Parameters of the virtual-bond valence angles:'
962 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
963 ' ATHETA0 ',' A1 ',' A2 ',&
966 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
967 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
969 write (iout,'(/a/9x,5a/79(1h-))') &
970 'Parameters of the expression for sigma(theta_c):',&
971 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
972 ' ALPH3 ',' SIGMA0C '
974 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
975 (polthet(j,i),j=0,3),sigc0(i)
977 write (iout,'(/a/9x,5a/79(1h-))') &
978 'Parameters of the second gaussian:',&
979 ' THETA0 ',' SIGMA0 ',' G1 ',&
982 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
983 sig0(i),(gthet(j,i),j=1,3)
987 'Parameters of the virtual-bond valence angles:'
988 write (iout,'(/a/9x,5a/79(1h-))') &
989 'Coefficients of expansion',&
990 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
991 ' b1*10^1 ',' b2*10^1 '
993 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
994 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
995 (10*bthet(j,i,1,1),j=1,2)
997 write (iout,'(/a/9x,5a/79(1h-))') &
998 'Parameters of the expression for sigma(theta_c):',&
999 ' alpha0 ',' alph1 ',' alph2 ',&
1000 ' alhp3 ',' sigma0c '
1002 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1003 (polthet(j,i),j=0,3),sigc0(i)
1005 write (iout,'(/a/9x,5a/79(1h-))') &
1006 'Parameters of the second gaussian:',&
1007 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1010 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1011 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1017 ! Read the parameters of Utheta determined from ab initio surfaces
1018 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1020 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1021 ntheterm3,nsingle,ndouble
1022 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1024 !----------------------------------------------------
1025 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1026 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1027 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1028 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1029 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1032 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1033 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1034 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1035 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1036 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1037 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1038 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1039 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1040 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1041 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1042 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1043 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1044 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1045 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1046 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1047 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1048 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1050 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1052 ithetyp(i)=-ithetyp(-i)
1055 aa0thet(:,:,:,:)=0.0d0
1056 aathet(:,:,:,:,:)=0.0d0
1057 bbthet(:,:,:,:,:,:)=0.0d0
1058 ccthet(:,:,:,:,:,:)=0.0d0
1059 ddthet(:,:,:,:,:,:)=0.0d0
1060 eethet(:,:,:,:,:,:)=0.0d0
1061 ffthet(:,:,:,:,:,:,:)=0.0d0
1062 ggthet(:,:,:,:,:,:,:)=0.0d0
1064 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1066 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1067 ! VAR:1=non-glicyne non-proline 2=proline
1068 ! VAR:negative values for D-aminoacid
1070 do j=-nthetyp,nthetyp
1071 do k=-nthetyp,nthetyp
1072 read (ithep,'(6a)',end=111,err=111) res1
1073 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1074 ! VAR: aa0thet is variable describing the average value of Foureir
1075 ! VAR: expansion series
1076 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1077 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1078 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1079 read (ithep,*,end=111,err=111) &
1080 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1081 read (ithep,*,end=111,err=111) &
1082 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1083 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1084 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1085 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1087 read (ithep,*,end=111,err=111) &
1088 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1089 ffthet(lll,llll,ll,i,j,k,iblock),&
1090 ggthet(llll,lll,ll,i,j,k,iblock),&
1091 ggthet(lll,llll,ll,i,j,k,iblock),&
1092 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1097 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1098 ! coefficients of theta-and-gamma-dependent terms are zero.
1099 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1100 ! RECOMENTDED AFTER VERSION 3.3)
1104 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1105 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1107 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1108 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1111 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1113 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1116 ! AND COMMENT THE LOOPS BELOW
1120 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1121 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1123 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1124 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1127 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1129 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1134 ! Substitution for D aminoacids from symmetry.
1137 do j=-nthetyp,nthetyp
1138 do k=-nthetyp,nthetyp
1139 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1141 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1145 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1146 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1147 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1148 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1154 ffthet(llll,lll,ll,i,j,k,iblock)= &
1155 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1156 ffthet(lll,llll,ll,i,j,k,iblock)= &
1157 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1158 ggthet(llll,lll,ll,i,j,k,iblock)= &
1159 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1160 ggthet(lll,llll,ll,i,j,k,iblock)= &
1161 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1170 ! Control printout of the coefficients of virtual-bond-angle potentials
1173 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1178 write (iout,'(//4a)') &
1179 'Type ',onelett(i),onelett(j),onelett(k)
1180 write (iout,'(//a,10x,a)') " l","a[l]"
1181 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1182 write (iout,'(i2,1pe15.5)') &
1183 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1185 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1186 "b",l,"c",l,"d",l,"e",l
1188 write (iout,'(i2,4(1pe15.5))') m,&
1189 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1190 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1194 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1195 "f+",l,"f-",l,"g+",l,"g-",l
1198 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1199 ffthet(n,m,l,i,j,k,iblock),&
1200 ffthet(m,n,l,i,j,k,iblock),&
1201 ggthet(n,m,l,i,j,k,iblock),&
1202 ggthet(m,n,l,i,j,k,iblock)
1212 write (2,*) "Start reading THETA_PDB",ithep_pdb
1214 ! write (2,*) 'i=',i
1215 read (ithep_pdb,*,err=111,end=111) &
1216 a0thet(i),(athet(j,i,1,1),j=1,2),&
1217 (bthet(j,i,1,1),j=1,2)
1218 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1219 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1220 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1221 sigc0(i)=sigc0(i)**2
1224 athet(1,i,1,-1)=athet(1,i,1,1)
1225 athet(2,i,1,-1)=athet(2,i,1,1)
1226 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1227 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1228 athet(1,i,-1,1)=-athet(1,i,1,1)
1229 athet(2,i,-1,1)=-athet(2,i,1,1)
1230 bthet(1,i,-1,1)=bthet(1,i,1,1)
1231 bthet(2,i,-1,1)=bthet(2,i,1,1)
1234 a0thet(i)=a0thet(-i)
1235 athet(1,i,-1,-1)=athet(1,-i,1,1)
1236 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1237 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1238 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1239 athet(1,i,-1,1)=athet(1,-i,1,1)
1240 athet(2,i,-1,1)=-athet(2,-i,1,1)
1241 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1242 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1243 athet(1,i,1,-1)=-athet(1,-i,1,1)
1244 athet(2,i,1,-1)=athet(2,-i,1,1)
1245 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1246 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1247 theta0(i)=theta0(-i)
1251 polthet(j,i)=polthet(j,-i)
1254 gthet(j,i)=gthet(j,-i)
1257 write (2,*) "End reading THETA_PDB"
1261 !--------------- Reading theta parameters for nucleic acid-------
1262 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1263 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1264 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1265 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1266 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1267 nthetyp_nucl+1,nthetyp_nucl+1))
1268 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1269 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1270 nthetyp_nucl+1,nthetyp_nucl+1))
1271 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1272 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1273 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1274 nthetyp_nucl+1,nthetyp_nucl+1))
1275 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1276 nthetyp_nucl+1,nthetyp_nucl+1))
1277 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1278 nthetyp_nucl+1,nthetyp_nucl+1))
1279 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1280 nthetyp_nucl+1,nthetyp_nucl+1))
1281 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1282 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1283 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1284 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1285 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1286 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1288 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1289 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1291 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1293 aa0thet_nucl(:,:,:)=0.0d0
1294 aathet_nucl(:,:,:,:)=0.0d0
1295 bbthet_nucl(:,:,:,:,:)=0.0d0
1296 ccthet_nucl(:,:,:,:,:)=0.0d0
1297 ddthet_nucl(:,:,:,:,:)=0.0d0
1298 eethet_nucl(:,:,:,:,:)=0.0d0
1299 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1300 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1305 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1306 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1307 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1308 read (ithep_nucl,*,end=111,err=111) &
1309 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1310 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1311 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1312 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1313 read (ithep_nucl,*,end=111,err=111) &
1314 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1315 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1316 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1321 !-------------------------------------------
1322 allocate(nlob(ntyp1)) !(ntyp1)
1323 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1324 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1325 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1332 gaussc(:,:,:,:)=0.0D0
1336 ! Read the parameters of the probability distribution/energy expression
1337 ! of the side chains.
1340 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1344 dsc_inv(i)=1.0D0/dsc(i)
1355 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1356 ((blower(k,l,1),l=1,k),k=1,3)
1357 censc(1,1,-i)=censc(1,1,i)
1358 censc(2,1,-i)=censc(2,1,i)
1359 censc(3,1,-i)=-censc(3,1,i)
1361 read (irotam,*,end=112,err=112) bsc(j,i)
1362 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1363 ((blower(k,l,j),l=1,k),k=1,3)
1364 censc(1,j,-i)=censc(1,j,i)
1365 censc(2,j,-i)=censc(2,j,i)
1366 censc(3,j,-i)=-censc(3,j,i)
1367 ! BSC is amplitude of Gaussian
1374 akl=akl+blower(k,m,j)*blower(l,m,j)
1378 if (((k.eq.3).and.(l.ne.3)) &
1379 .or.((l.eq.3).and.(k.ne.3))) then
1380 gaussc(k,l,j,-i)=-akl
1381 gaussc(l,k,j,-i)=-akl
1383 gaussc(k,l,j,-i)=akl
1384 gaussc(l,k,j,-i)=akl
1393 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1396 if (nlobi.gt.0) then
1398 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1399 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1400 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1401 'log h',(bsc(j,i),j=1,nlobi)
1402 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1403 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1405 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1406 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1409 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1410 write (iout,'(a,f10.4,4(16x,f10.4))') &
1411 'Center ',(bsc(j,i),j=1,nlobi)
1412 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1421 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1422 ! added by Urszula Kozlowska 07/11/2007
1424 !el Maximum number of SC local term fitting function coefficiants
1425 !el integer,parameter :: maxsccoef=65
1427 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1430 read (irotam,*,end=112,err=112)
1432 read (irotam,*,end=112,err=112)
1435 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1439 !---------reading nucleic acid parameters for rotamers-------------------
1440 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1441 do i=1,ntyp_molec(2)
1442 read (irotam_nucl,*,end=112,err=112)
1444 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1450 write (iout,*) "Base rotamer parameters"
1451 do i=1,ntyp_molec(2)
1452 write (iout,'(a)') restyp(i,2)
1453 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1458 ! Read the parameters of the probability distribution/energy expression
1459 ! of the side chains.
1461 write (2,*) "Start reading ROTAM_PDB"
1463 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1467 dsc_inv(i)=1.0D0/dsc(i)
1478 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1479 ((blower(k,l,1),l=1,k),k=1,3)
1481 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1482 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1483 ((blower(k,l,j),l=1,k),k=1,3)
1490 akl=akl+blower(k,m,j)*blower(l,m,j)
1500 write (2,*) "End reading ROTAM_PDB"
1506 ! Read torsional parameters in old format
1508 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1510 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1511 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1512 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1514 !el from energy module--------
1515 allocate(v1(nterm_old,ntortyp,ntortyp))
1516 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1517 !el---------------------------
1522 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1528 write (iout,'(/a/)') 'Torsional constants:'
1531 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1532 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1538 ! Read torsional parameters
1540 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1541 read (itorp,*,end=113,err=113) ntortyp
1542 !el from energy module---------
1543 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1544 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1546 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1547 allocate(vlor2(maxlor,ntortyp,ntortyp))
1548 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1549 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1551 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1552 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1553 !el---------------------------
1556 !el---------------------------
1558 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1560 itortyp(i)=-itortyp(-i)
1562 itortyp(ntyp1)=ntortyp
1563 itortyp(-ntyp1)=-ntortyp
1565 write (iout,*) 'ntortyp',ntortyp
1567 do j=-ntortyp+1,ntortyp-1
1568 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1570 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1571 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1574 do k=1,nterm(i,j,iblock)
1575 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1577 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1578 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1579 v0ij=v0ij+si*v1(k,i,j,iblock)
1581 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1582 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1583 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1585 do k=1,nlor(i,j,iblock)
1586 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1587 vlor2(k,i,j),vlor3(k,i,j)
1588 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1591 v0(-i,-j,iblock)=v0ij
1597 write (iout,'(/a/)') 'Torsional constants:'
1599 do i=-ntortyp,ntortyp
1600 do j=-ntortyp,ntortyp
1601 write (iout,*) 'ityp',i,' jtyp',j
1602 write (iout,*) 'Fourier constants'
1603 do k=1,nterm(i,j,iblock)
1604 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1607 write (iout,*) 'Lorenz constants'
1608 do k=1,nlor(i,j,iblock)
1609 write (iout,'(3(1pe15.5))') &
1610 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1616 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1618 ! 6/23/01 Read parameters for double torsionals
1620 !el from energy module------------
1621 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1622 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1623 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1624 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1625 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1626 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1627 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1628 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1629 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1630 !---------------------------------
1634 do j=-ntortyp+1,ntortyp-1
1635 do k=-ntortyp+1,ntortyp-1
1636 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1637 ! write (iout,*) "OK onelett",
1640 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1641 .or. t3.ne.toronelet(k)) then
1642 write (iout,*) "Error in double torsional parameter file",&
1645 call MPI_Finalize(Ierror)
1647 stop "Error in double torsional parameter file"
1649 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1650 ntermd_2(i,j,k,iblock)
1651 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1652 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1653 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1654 ntermd_1(i,j,k,iblock))
1655 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1656 ntermd_1(i,j,k,iblock))
1657 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1658 ntermd_1(i,j,k,iblock))
1659 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1660 ntermd_1(i,j,k,iblock))
1661 ! Martix of D parameters for one dimesional foureir series
1662 do l=1,ntermd_1(i,j,k,iblock)
1663 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1664 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1665 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1666 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1667 ! write(iout,*) "whcodze" ,
1668 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1670 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1671 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1672 v2s(m,l,i,j,k,iblock),&
1673 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1674 ! Martix of D parameters for two dimesional fourier series
1675 do l=1,ntermd_2(i,j,k,iblock)
1677 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1678 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1679 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1680 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1689 write (iout,*) 'Constants for double torsionals'
1692 do j=-ntortyp+1,ntortyp-1
1693 do k=-ntortyp+1,ntortyp-1
1694 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1695 ' nsingle',ntermd_1(i,j,k,iblock),&
1696 ' ndouble',ntermd_2(i,j,k,iblock)
1698 write (iout,*) 'Single angles:'
1699 do l=1,ntermd_1(i,j,k,iblock)
1700 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1701 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1702 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1703 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1706 write (iout,*) 'Pairs of angles:'
1707 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1708 do l=1,ntermd_2(i,j,k,iblock)
1709 write (iout,'(i5,20f10.5)') &
1710 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1713 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1714 do l=1,ntermd_2(i,j,k,iblock)
1715 write (iout,'(i5,20f10.5)') &
1716 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1717 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1726 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1727 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1728 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1729 !el from energy module---------
1730 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1731 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1733 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1734 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1735 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1736 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1738 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1739 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1740 !el---------------------------
1743 !el--------------------
1744 read (itorp_nucl,*,end=113,err=113) &
1745 (itortyp_nucl(i),i=1,ntyp_molec(2))
1746 ! print *,itortyp_nucl(:)
1747 !c write (iout,*) 'ntortyp',ntortyp
1750 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1751 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1754 do k=1,nterm_nucl(i,j)
1755 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1756 v0ij=v0ij+si*v1_nucl(k,i,j)
1759 do k=1,nlor_nucl(i,j)
1760 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1761 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1762 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1768 ! Read of Side-chain backbone correlation parameters
1769 ! Modified 11 May 2012 by Adasko
1772 read (isccor,*,end=119,err=119) nsccortyp
1774 !el from module energy-------------
1775 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1776 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1777 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1778 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1779 !-----------------------------------
1781 !el from module energy-------------
1782 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1784 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1786 isccortyp(i)=-isccortyp(-i)
1788 iscprol=isccortyp(20)
1789 ! write (iout,*) 'ntortyp',ntortyp
1791 !c maxinter is maximum interaction sites
1792 !el from module energy---------
1793 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1794 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1795 -nsccortyp:nsccortyp))
1796 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1797 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1798 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1799 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1800 !-----------------------------------
1804 read (isccor,*,end=119,err=119) &
1805 nterm_sccor(i,j),nlor_sccor(i,j)
1811 nterm_sccor(-i,j)=nterm_sccor(i,j)
1812 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1813 nterm_sccor(i,-j)=nterm_sccor(i,j)
1814 do k=1,nterm_sccor(i,j)
1815 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1817 if (j.eq.iscprol) then
1818 if (i.eq.isccortyp(10)) then
1819 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1820 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1822 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1823 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1824 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1825 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1826 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1827 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1828 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1829 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1832 if (i.eq.isccortyp(10)) then
1833 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1834 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1836 if (j.eq.isccortyp(10)) then
1837 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1838 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1840 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1841 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1842 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1843 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1844 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1845 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1849 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1850 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1851 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1852 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1855 do k=1,nlor_sccor(i,j)
1856 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1857 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1858 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1859 (1+vlor3sccor(k,i,j)**2)
1861 v0sccor(l,i,j)=v0ijsccor
1862 v0sccor(l,-i,j)=v0ijsccor1
1863 v0sccor(l,i,-j)=v0ijsccor2
1864 v0sccor(l,-i,-j)=v0ijsccor3
1870 !el from module energy-------------
1871 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1873 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1874 ! write (iout,*) 'ntortyp',ntortyp
1876 !c maxinter is maximum interaction sites
1877 !el from module energy---------
1878 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1879 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1880 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1881 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1882 !-----------------------------------
1886 read (isccor,*,end=119,err=119) &
1887 nterm_sccor(i,j),nlor_sccor(i,j)
1891 do k=1,nterm_sccor(i,j)
1892 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1894 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1897 do k=1,nlor_sccor(i,j)
1898 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1899 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1900 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1901 (1+vlor3sccor(k,i,j)**2)
1903 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1912 write (iout,'(/a/)') 'Torsional constants:'
1915 write (iout,*) 'ityp',i,' jtyp',j
1916 write (iout,*) 'Fourier constants'
1917 do k=1,nterm_sccor(i,j)
1918 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1920 write (iout,*) 'Lorenz constants'
1921 do k=1,nlor_sccor(i,j)
1922 write (iout,'(3(1pe15.5))') &
1923 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1930 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1931 ! interaction energy of the Gly, Ala, and Pro prototypes.
1936 write (iout,*) "Coefficients of the cumulants"
1938 read (ifourier,*) nloctyp
1939 !write(iout,*) "nloctyp",nloctyp
1940 !el from module energy-------
1941 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1942 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1943 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1944 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1945 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1946 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1947 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1948 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1952 !--------------------------------
1955 read (ifourier,*,end=115,err=115)
1956 read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
1958 write (iout,*) 'Type',i
1959 write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
1967 B1tilde(1,i) = buse(3)
1968 B1tilde(2,i) =-buse(5)
1969 B1tilde(1,-i) =-buse(3)
1970 B1tilde(2,-i) =buse(5)
1971 ! buse1tilde(1,i)=0.0d0
1972 ! buse1tilde(2,i)=0.0d0
1992 Ctilde(1,1,i)=buse(7)
1993 Ctilde(1,2,i)=buse(9)
1994 Ctilde(2,1,i)=-buse(9)
1995 Ctilde(2,2,i)=buse(7)
1996 Ctilde(1,1,-i)=buse(7)
1997 Ctilde(1,2,-i)=-buse(9)
1998 Ctilde(2,1,-i)=buse(9)
1999 Ctilde(2,2,-i)=buse(7)
2001 ! Ctilde(1,1,i)=0.0d0
2002 ! Ctilde(1,2,i)=0.0d0
2003 ! Ctilde(2,1,i)=0.0d0
2004 ! Ctilde(2,2,i)=0.0d0
2017 Dtilde(1,1,i)=buse(6)
2018 Dtilde(1,2,i)=buse(8)
2019 Dtilde(2,1,i)=-buse(8)
2020 Dtilde(2,2,i)=buse(6)
2021 Dtilde(1,1,-i)=buse(6)
2022 Dtilde(1,2,-i)=-buse(8)
2023 Dtilde(2,1,-i)=buse(8)
2024 Dtilde(2,2,-i)=buse(6)
2026 ! Dtilde(1,1,i)=0.0d0
2027 ! Dtilde(1,2,i)=0.0d0
2028 ! Dtilde(2,1,i)=0.0d0
2029 ! Dtilde(2,2,i)=0.0d0
2030 EE(1,1,i)= buse(10)+buse(11)
2031 EE(2,2,i)=-buse(10)+buse(11)
2032 EE(2,1,i)= buse(12)-buse(13)
2033 EE(1,2,i)= buse(12)+buse(13)
2034 EE(1,1,-i)= buse(10)+buse(11)
2035 EE(2,2,-i)=-buse(10)+buse(11)
2036 EE(2,1,-i)=-buse(12)+buse(13)
2037 EE(1,2,-i)=-buse(12)-buse(13)
2043 ! ee(2,1,i)=ee(1,2,i)
2047 write (iout,*) 'Type',i
2049 write(iout,*) B1(1,i),B1(2,i)
2051 write(iout,*) B2(1,i),B2(2,i)
2054 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2058 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2062 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2067 ! Read electrostatic-interaction parameters
2072 write (iout,'(/a)') 'Electrostatic interaction constants:'
2073 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2074 'IT','JT','APP','BPP','AEL6','AEL3'
2076 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2077 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2078 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2079 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2084 app (i,j)=epp(i,j)*rri*rri
2085 bpp (i,j)=-2.0D0*epp(i,j)*rri
2086 ael6(i,j)=elpp6(i,j)*4.2D0**6
2087 ael3(i,j)=elpp3(i,j)*4.2D0**3
2089 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2095 ! Read side-chain interaction parameters.
2097 !el from module energy - COMMON.INTERACT-------
2098 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2099 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2100 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2102 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2103 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2104 allocate(epslip(ntyp,ntyp))
2112 !--------------------------------
2114 read (isidep,*,end=117,err=117) ipot,expon
2115 if (ipot.lt.1 .or. ipot.gt.5) then
2116 write (iout,'(2a)') 'Error while reading SC interaction',&
2117 'potential file - unknown potential type.'
2119 call MPI_Finalize(Ierror)
2125 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2126 ', exponents are ',expon,2*expon
2127 ! goto (10,20,30,30,40) ipot
2129 !----------------------- LJ potential ---------------------------------
2131 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2132 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2133 (sigma0(i),i=1,ntyp)
2135 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2136 write (iout,'(a/)') 'The epsilon array:'
2137 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2138 write (iout,'(/a)') 'One-body parameters:'
2139 write (iout,'(a,4x,a)') 'residue','sigma'
2140 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2143 !----------------------- LJK potential --------------------------------
2145 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2146 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2147 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2149 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2150 write (iout,'(a/)') 'The epsilon array:'
2151 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2152 write (iout,'(/a)') 'One-body parameters:'
2153 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2154 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2158 !---------------------- GB or BP potential -----------------------------
2161 ! print *,"I AM in SCELE",scelemode
2162 if (scelemode.eq.0) then
2164 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2166 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2167 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2168 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2169 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2171 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2174 ! For the GB potential convert sigma'**2 into chi'
2177 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2181 write (iout,'(/a/)') 'Parameters of the BP potential:'
2182 write (iout,'(a/)') 'The epsilon array:'
2183 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2184 write (iout,'(/a)') 'One-body parameters:'
2185 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2187 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2188 chip(i),alp(i),i=1,ntyp)
2191 ! print *,ntyp,"NTYP"
2192 allocate(icharge(ntyp1))
2193 ! print *,ntyp,icharge(i)
2195 read (isidep,*) (icharge(i),i=1,ntyp)
2196 print *,ntyp,icharge(i)
2197 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2198 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2199 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2200 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2201 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2202 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2203 allocate(epsintab(ntyp,ntyp))
2204 allocate(dtail(2,ntyp,ntyp))
2205 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2206 allocate(wqdip(2,ntyp,ntyp))
2207 allocate(wstate(4,ntyp,ntyp))
2208 allocate(dhead(2,2,ntyp,ntyp))
2209 allocate(nstate(ntyp,ntyp))
2210 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2211 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2214 ! write (*,*) "Im in ALAB", i, " ", j
2216 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2217 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2218 chis(i,j),chis(j,i), &
2219 nstate(i,j),(wstate(k,i,j),k=1,4), &
2220 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2221 dtail(1,i,j),dtail(2,i,j), &
2222 epshead(i,j),sig0head(i,j), &
2223 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2224 alphapol(i,j),alphapol(j,i), &
2225 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
2226 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2232 sigma(i,j) = sigma(j,i)
2233 sigmap1(i,j)=sigmap1(j,i)
2234 sigmap2(i,j)=sigmap2(j,i)
2235 sigiso1(i,j)=sigiso1(j,i)
2236 sigiso2(i,j)=sigiso2(j,i)
2237 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2238 nstate(i,j) = nstate(j,i)
2239 dtail(1,i,j) = dtail(1,j,i)
2240 dtail(2,i,j) = dtail(2,j,i)
2242 alphasur(k,i,j) = alphasur(k,j,i)
2243 wstate(k,i,j) = wstate(k,j,i)
2244 alphiso(k,i,j) = alphiso(k,j,i)
2247 dhead(2,1,i,j) = dhead(1,1,j,i)
2248 dhead(2,2,i,j) = dhead(1,2,j,i)
2249 dhead(1,1,i,j) = dhead(2,1,j,i)
2250 dhead(1,2,i,j) = dhead(2,2,j,i)
2252 epshead(i,j) = epshead(j,i)
2253 sig0head(i,j) = sig0head(j,i)
2256 wqdip(k,i,j) = wqdip(k,j,i)
2259 wquad(i,j) = wquad(j,i)
2260 epsintab(i,j) = epsintab(j,i)
2261 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2266 !--------------------- GBV potential -----------------------------------
2268 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2269 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2270 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2271 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2273 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2274 write (iout,'(a/)') 'The epsilon array:'
2275 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2276 write (iout,'(/a)') 'One-body parameters:'
2277 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2278 's||/s_|_^2',' chip ',' alph '
2279 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2280 sigii(i),chip(i),alp(i),i=1,ntyp)
2283 write(iout,*)"Wrong ipot"
2289 !-----------------------------------------------------------------------
2290 ! Calculate the "working" parameters of SC interactions.
2292 !el from module energy - COMMON.INTERACT-------
2293 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2294 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2295 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2296 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2297 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2298 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2300 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2306 if (scelemode.eq.0) then
2315 sc_aa_tube_par(:)=0.0d0
2316 sc_bb_tube_par(:)=0.0d0
2318 !--------------------------------
2323 epslip(i,j)=epslip(j,i)
2326 if (scelemode.eq.0) then
2329 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2330 sigma(j,i)=sigma(i,j)
2331 rs0(i,j)=dwa16*sigma(i,j)
2336 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2337 'Working parameters of the SC interactions:',&
2338 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2343 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2345 ! print *,"SIGMA ZLA?",sigma(i,j)
2353 sigeps=dsign(1.0D0,epsij)
2355 aa_aq(i,j)=epsij*rrij*rrij
2356 print *,"ADASKO",epsij,rrij,expon
2357 bb_aq(i,j)=-sigeps*epsij*rrij
2358 aa_aq(j,i)=aa_aq(i,j)
2359 bb_aq(j,i)=bb_aq(i,j)
2360 epsijlip=epslip(i,j)
2361 sigeps=dsign(1.0D0,epsijlip)
2362 epsijlip=dabs(epsijlip)
2363 aa_lip(i,j)=epsijlip*rrij*rrij
2364 bb_lip(i,j)=-sigeps*epsijlip*rrij
2365 aa_lip(j,i)=aa_lip(i,j)
2366 bb_lip(j,i)=bb_lip(i,j)
2367 !C write(iout,*) aa_lip
2368 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2369 sigt1sq=sigma0(i)**2
2370 sigt2sq=sigma0(j)**2
2373 ratsig1=sigt2sq/sigt1sq
2374 ratsig2=1.0D0/ratsig1
2375 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2376 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2377 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2381 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2382 sigmaii(i,j)=rsum_max
2383 sigmaii(j,i)=rsum_max
2385 ! sigmaii(i,j)=r0(i,j)
2386 ! sigmaii(j,i)=r0(i,j)
2388 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2389 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2390 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2391 augm(i,j)=epsij*r_augm**(2*expon)
2392 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2399 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2400 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2401 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2406 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2407 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2408 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2409 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2410 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2411 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2412 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2413 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2414 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2415 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2416 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2417 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2418 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2419 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2420 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2421 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2430 read (isidep_nucl,*) ipot_nucl
2431 ! print *,"TU?!",ipot_nucl
2432 if (ipot_nucl.eq.1) then
2433 do i=1,ntyp_molec(2)
2434 do j=i,ntyp_molec(2)
2435 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2436 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2440 do i=1,ntyp_molec(2)
2441 do j=i,ntyp_molec(2)
2442 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2443 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2444 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2448 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2449 do i=1,ntyp_molec(2)
2450 do j=i,ntyp_molec(2)
2451 rrij=sigma_nucl(i,j)
2455 epsij=4*eps_nucl(i,j)
2456 sigeps=dsign(1.0D0,epsij)
2458 aa_nucl(i,j)=epsij*rrij*rrij
2459 bb_nucl(i,j)=-sigeps*epsij*rrij
2460 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2461 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2462 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2463 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2464 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2465 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2468 aa_nucl(i,j)=aa_nucl(j,i)
2469 bb_nucl(i,j)=bb_nucl(j,i)
2470 ael3_nucl(i,j)=ael3_nucl(j,i)
2471 ael6_nucl(i,j)=ael6_nucl(j,i)
2472 ael63_nucl(i,j)=ael63_nucl(j,i)
2473 ael32_nucl(i,j)=ael32_nucl(j,i)
2474 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2475 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2476 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2477 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2478 eps_nucl(i,j)=eps_nucl(j,i)
2479 sigma_nucl(i,j)=sigma_nucl(j,i)
2480 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2484 write(iout,*) "tube param"
2485 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2486 ccavtubpep,dcavtubpep,tubetranenepep
2487 sigmapeptube=sigmapeptube**6
2488 sigeps=dsign(1.0D0,epspeptube)
2489 epspeptube=dabs(epspeptube)
2490 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2491 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2492 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2494 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2495 ccavtub(i),dcavtub(i),tubetranene(i)
2496 sigmasctube=sigmasctube**6
2497 sigeps=dsign(1.0D0,epssctube)
2498 epssctube=dabs(epssctube)
2499 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2500 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2501 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2503 !-----------------READING SC BASE POTENTIALS-----------------------------
2504 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2505 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2506 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2507 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2508 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2509 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2510 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2511 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2512 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2513 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2514 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2515 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2516 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2517 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2518 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2519 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2522 do i=1,ntyp_molec(1)
2523 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2524 write (*,*) "Im in ", i, " ", j
2525 read(isidep_scbase,*) &
2526 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2527 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2528 write(*,*) "eps",eps_scbase(i,j)
2529 read(isidep_scbase,*) &
2530 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2531 chis_scbase(i,j,1),chis_scbase(i,j,2)
2532 read(isidep_scbase,*) &
2533 dhead_scbasei(i,j), &
2534 dhead_scbasej(i,j), &
2535 rborn_scbasei(i,j),rborn_scbasej(i,j)
2536 read(isidep_scbase,*) &
2537 (wdipdip_scbase(k,i,j),k=1,3), &
2538 (wqdip_scbase(k,i,j),k=1,2)
2539 read(isidep_scbase,*) &
2540 alphapol_scbase(i,j), &
2541 epsintab_scbase(i,j)
2544 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2545 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2547 do i=1,ntyp_molec(1)
2548 do j=1,ntyp_molec(2)-1
2549 epsij=eps_scbase(i,j)
2550 rrij=sigma_scbase(i,j)
2555 sigeps=dsign(1.0D0,epsij)
2557 aa_scbase(i,j)=epsij*rrij*rrij
2558 bb_scbase(i,j)=-sigeps*epsij*rrij
2561 !-----------------READING PEP BASE POTENTIALS-------------------
2562 allocate(eps_pepbase(ntyp_molec(2)))
2563 allocate(sigma_pepbase(ntyp_molec(2)))
2564 allocate(chi_pepbase(ntyp_molec(2),2))
2565 allocate(chipp_pepbase(ntyp_molec(2),2))
2566 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2567 allocate(sigmap1_pepbase(ntyp_molec(2)))
2568 allocate(sigmap2_pepbase(ntyp_molec(2)))
2569 allocate(chis_pepbase(ntyp_molec(2),2))
2570 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2573 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2574 write (*,*) "Im in ", i, " ", j
2575 read(isidep_pepbase,*) &
2576 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2577 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2578 write(*,*) "eps",eps_pepbase(j)
2579 read(isidep_pepbase,*) &
2580 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2581 chis_pepbase(j,1),chis_pepbase(j,2)
2582 read(isidep_pepbase,*) &
2583 (wdipdip_pepbase(k,j),k=1,3)
2585 allocate(aa_pepbase(ntyp_molec(2)))
2586 allocate(bb_pepbase(ntyp_molec(2)))
2588 do j=1,ntyp_molec(2)-1
2589 epsij=eps_pepbase(j)
2590 rrij=sigma_pepbase(j)
2595 sigeps=dsign(1.0D0,epsij)
2597 aa_pepbase(j)=epsij*rrij*rrij
2598 bb_pepbase(j)=-sigeps*epsij*rrij
2600 !--------------READING SC PHOSPHATE-------------------------------------
2601 allocate(eps_scpho(ntyp_molec(1)))
2602 allocate(sigma_scpho(ntyp_molec(1)))
2603 allocate(chi_scpho(ntyp_molec(1),2))
2604 allocate(chipp_scpho(ntyp_molec(1),2))
2605 allocate(alphasur_scpho(4,ntyp_molec(1)))
2606 allocate(sigmap1_scpho(ntyp_molec(1)))
2607 allocate(sigmap2_scpho(ntyp_molec(1)))
2608 allocate(chis_scpho(ntyp_molec(1),2))
2609 allocate(wqq_scpho(ntyp_molec(1)))
2610 allocate(wqdip_scpho(2,ntyp_molec(1)))
2611 allocate(alphapol_scpho(ntyp_molec(1)))
2612 allocate(epsintab_scpho(ntyp_molec(1)))
2613 allocate(dhead_scphoi(ntyp_molec(1)))
2614 allocate(rborn_scphoi(ntyp_molec(1)))
2615 allocate(rborn_scphoj(ntyp_molec(1)))
2616 allocate(alphi_scpho(ntyp_molec(1)))
2620 do j=1,ntyp_molec(1) ! without U then we will take T for U
2621 write (*,*) "Im in scpho ", i, " ", j
2622 read(isidep_scpho,*) &
2623 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2624 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2625 write(*,*) "eps",eps_scpho(j)
2626 read(isidep_scpho,*) &
2627 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2628 chis_scpho(j,1),chis_scpho(j,2)
2629 read(isidep_scpho,*) &
2630 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2631 read(isidep_scpho,*) &
2632 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2636 allocate(aa_scpho(ntyp_molec(1)))
2637 allocate(bb_scpho(ntyp_molec(1)))
2639 do j=1,ntyp_molec(1)
2646 sigeps=dsign(1.0D0,epsij)
2648 aa_scpho(j)=epsij*rrij*rrij
2649 bb_scpho(j)=-sigeps*epsij*rrij
2653 read(isidep_peppho,*) &
2654 eps_peppho,sigma_peppho
2655 read(isidep_peppho,*) &
2656 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2657 read(isidep_peppho,*) &
2658 (wqdip_peppho(k),k=1,2)
2666 sigeps=dsign(1.0D0,epsij)
2668 aa_peppho=epsij*rrij*rrij
2669 bb_peppho=-sigeps*epsij*rrij
2672 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2677 ! Define the SC-p interaction constants (hard-coded; old style)
2680 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2682 ! aad(i,1)=0.3D0*4.0D0**12
2683 ! Following line for constants currently implemented
2684 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2685 aad(i,1)=1.5D0*4.0D0**12
2686 ! aad(i,1)=0.17D0*5.6D0**12
2688 ! "Soft" SC-p repulsion
2690 ! Following line for constants currently implemented
2691 ! aad(i,1)=0.3D0*4.0D0**6
2692 ! "Hard" SC-p repulsion
2693 bad(i,1)=3.0D0*4.0D0**6
2694 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2703 ! 8/9/01 Read the SC-p interaction constants from file
2706 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2709 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2710 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2711 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2712 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2716 write (iout,*) "Parameters of SC-p interactions:"
2718 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2719 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2724 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2726 do i=1,ntyp_molec(2)
2727 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2729 do i=1,ntyp_molec(2)
2730 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2731 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2733 r0pp=1.12246204830937298142*5.16158
2739 ! Define the constants of the disulfide bridge
2743 ! Old arbitrary potential - commented out.
2748 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2749 ! energy surface of diethyl disulfide.
2750 ! A. Liwo and U. Kozlowska, 11/24/03
2767 write (iout,'(/a)') "Disulfide bridge parameters:"
2768 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2769 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2770 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2771 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2775 111 write (iout,*) "Error reading bending energy parameters."
2777 112 write (iout,*) "Error reading rotamer energy parameters."
2779 113 write (iout,*) "Error reading torsional energy parameters."
2781 114 write (iout,*) "Error reading double torsional energy parameters."
2783 115 write (iout,*) &
2784 "Error reading cumulant (multibody energy) parameters."
2786 116 write (iout,*) "Error reading electrostatic energy parameters."
2788 117 write (iout,*) "Error reading side chain interaction parameters."
2790 118 write (iout,*) "Error reading SCp interaction parameters."
2792 119 write (iout,*) "Error reading SCCOR parameters"
2795 call MPI_Finalize(Ierror)
2799 end subroutine parmread
2801 !-----------------------------------------------------------------------------
2803 !-----------------------------------------------------------------------------
2804 subroutine printmat(ldim,m,n,iout,key,a)
2807 character(len=3),dimension(n) :: key
2808 real(kind=8),dimension(ldim,n) :: a
2810 integer :: i,j,k,m,iout,nlim
2814 write (iout,1000) (key(k),k=i,nlim)
2816 1000 format (/5x,8(6x,a3))
2817 1020 format (/80(1h-)/)
2819 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2822 1010 format (a3,2x,8(f9.4))
2824 end subroutine printmat
2825 !-----------------------------------------------------------------------------
2827 !-----------------------------------------------------------------------------
2829 ! Read the PDB file and convert the peptide geometry into virtual-chain
2832 use energy_data, only: itype,istype
2836 ! use control, only: rescode,sugarcode
2837 ! implicit real*8 (a-h,o-z)
2838 ! include 'DIMENSIONS'
2839 ! include 'COMMON.LOCAL'
2840 ! include 'COMMON.VAR'
2841 ! include 'COMMON.CHAIN'
2842 ! include 'COMMON.INTERACT'
2843 ! include 'COMMON.IOUNITS'
2844 ! include 'COMMON.GEO'
2845 ! include 'COMMON.NAMES'
2846 ! include 'COMMON.CONTROL'
2847 ! include 'COMMON.DISTFIT'
2848 ! include 'COMMON.SETUP'
2849 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2851 logical :: lprn=.true.,fail
2852 real(kind=8),dimension(3) :: e1,e2,e3
2853 real(kind=8) :: dcj,efree_temp
2854 character(len=3) :: seq,res
2855 character(len=5) :: atom
2856 character(len=80) :: card
2857 real(kind=8),dimension(3,20) :: sccor
2858 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2859 integer :: isugar,molecprev,firstion
2860 character*1 :: sugar
2862 real(kind=8),dimension(3) :: ccc
2864 integer,dimension(2,maxres/3) :: hfrag_alloc
2865 integer,dimension(4,maxres/3) :: bfrag_alloc
2866 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2867 real(kind=8),dimension(:,:), allocatable :: c_temporary
2868 integer,dimension(:,:) , allocatable :: itype_temporary
2869 integer,dimension(:),allocatable :: istype_temp
2876 ! write (2,*) "UNRES_PDB",unres_pdb
2896 !-----------------------------
2897 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2898 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2899 if(.not. allocated(istype)) allocate(istype(maxres))
2901 read (ipdbin,'(a80)',end=10) card
2902 ! write (iout,'(a)') card
2903 if (card(:5).eq.'HELIX') then
2906 read(card(22:25),*) hfrag(1,nhfrag)
2907 read(card(34:37),*) hfrag(2,nhfrag)
2909 if (card(:5).eq.'SHEET') then
2912 read(card(24:26),*) bfrag(1,nbfrag)
2913 read(card(35:37),*) bfrag(2,nbfrag)
2914 !rc----------------------------------------
2915 !rc to be corrected !!!
2916 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2917 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2918 !rc----------------------------------------
2920 if (card(:3).eq.'END') then
2922 else if (card(:3).eq.'TER') then
2927 itype(ires_old,molecule)=ntyp1_molec(molecule)
2928 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2929 nres_molec(molecule)=nres_molec(molecule)+2
2931 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2934 dc(j,ires)=sccor(j,iii)
2937 call sccenter(ires,iii,sccor)
2943 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2944 ! Fish out the ATOM cards.
2945 ! write(iout,*) 'card',card(1:20)
2946 if (index(card(1:4),'ATOM').gt.0) then
2947 read (card(12:16),*) atom
2948 ! write (iout,*) "! ",atom," !",ires
2949 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2950 read (card(23:26),*) ires
2951 read (card(18:20),'(a3)') res
2952 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2953 ! & " ires_old",ires_old
2954 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2955 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2956 if (ires-ishift+ishift1.ne.ires_old) then
2957 ! Calculate the CM of the preceding residue.
2958 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2960 ! write (iout,*) "Calculating sidechain center iii",iii
2963 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2966 call sccenter(ires_old,iii,sccor)
2970 ! Start new residue.
2971 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2974 else if (ibeg.eq.1) then
2975 write (iout,*) "BEG ires",ires
2977 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2980 nres_molec(molecule)=nres_molec(molecule)+1
2982 ires=ires-ishift+ishift1
2984 ! write (iout,*) "ishift",ishift," ires",ires,&
2985 ! " ires_old",ires_old
2987 else if (ibeg.eq.2) then
2989 ishift=-ires_old+ires-1 !!!!!
2990 ishift1=ishift1-1 !!!!!
2991 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2992 ires=ires-ishift+ishift1
2993 ! print *,ires,ishift,ishift1
2997 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2998 ires=ires-ishift+ishift1
3001 ! print *,'atom',ires,atom
3002 if (res.eq.'ACE' .or. res.eq.'NHE') then
3005 if (atom.eq.'CA '.or.atom.eq.'N ') then
3007 itype(ires,molecule)=rescode(ires,res,0,molecule)
3009 ! nres_molec(molecule)=nres_molec(molecule)+1
3012 itype(ires,molecule)=rescode(ires,res(2:3),0,molecule)
3013 ! nres_molec(molecule)=nres_molec(molecule)+1
3014 read (card(19:19),'(a1)') sugar
3015 isugar=sugarcode(sugar,ires)
3016 ! if (ibeg.eq.1) then
3020 ! print *,"ires=",ires,istype(ires)
3026 ires=ires-ishift+ishift1
3028 ! write (iout,*) "ires_old",ires_old," ires",ires
3029 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3032 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3033 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3034 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3035 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3036 ! print *,ires,ishift,ishift1
3037 ! write (iout,*) "backbone ",atom
3039 write (iout,'(2i3,2x,a,3f8.3)') &
3040 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3043 nres_molec(molecule)=nres_molec(molecule)+1
3045 sccor(j,iii)=c(j,ires)
3047 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3048 atom.eq."C2'" .or. atom.eq."C3'" &
3049 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3050 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3051 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3052 ! print *,ires,ishift,ishift1
3056 ! sccor(j,iii)=c(j,ires)
3059 c(j,ires)=c(j,ires)+ccc(j)/5.0
3061 print *,counter,molecule
3062 if (counter.eq.5) then
3064 nres_molec(molecule)=nres_molec(molecule)+1
3067 ! sccor(j,iii)=c(j,ires)
3071 ! print *, "ATOM",atom(1:3)
3072 else if (atom.eq."C5'") then
3073 read (card(19:19),'(a1)') sugar
3074 isugar=sugarcode(sugar,ires)
3079 ! print *,ires,istype(ires)
3083 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3084 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3085 nres_molec(molecule)=nres_molec(molecule)+1
3086 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3090 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3092 else if ((atom.eq."C1'").and.unres_pdb) then
3094 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3095 ! write (*,*) card(23:27),ires,itype(ires,1)
3096 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3097 atom.ne.'N' .and. atom.ne.'C' .and. &
3098 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3099 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3100 .and. atom.ne.'P '.and. &
3101 atom(1:1).ne.'H' .and. &
3102 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3104 ! write (iout,*) "sidechain ",atom
3105 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3106 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3107 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3109 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3112 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3113 if (firstion.eq.0) then
3117 dc(j,ires)=sccor(j,iii)
3120 call sccenter(ires,iii,sccor)
3123 read (card(12:16),*) atom
3124 print *,"HETATOM", atom
3125 read (card(18:20),'(a3)') res
3126 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3127 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3128 .or.(atom(1:2).eq.'K ')) &
3131 if (molecule.ne.5) molecprev=molecule
3133 nres_molec(molecule)=nres_molec(molecule)+1
3134 print *,"HERE",nres_molec(molecule)
3136 itype(ires,molecule)=rescode(ires,res,0,molecule)
3137 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3141 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3142 if (ires.eq.0) return
3143 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3146 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3148 nres_molec(molecule)=nres_molec(molecule)-2
3149 print *,'I have',nres, nres_molec(:)
3151 do k=1,4 ! ions are without dummy
3152 if (nres_molec(k).eq.0) cycle
3154 ! write (iout,*) i,itype(i,1)
3155 ! if (itype(i,1).eq.ntyp1) then
3156 ! write (iout,*) "dummy",i,itype(i,1)
3158 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3159 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3163 if (itype(i,k).eq.ntyp1_molec(k)) then
3164 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3165 if (itype(i+2,k).eq.0) then
3166 ! print *,"masz sieczke"
3168 if (itype(i+2,j).ne.0) then
3170 itype(i+1,j)=ntyp1_molec(j)
3171 nres_molec(k)=nres_molec(k)-1
3172 nres_molec(j)=nres_molec(j)+1
3178 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3179 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3180 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3181 ! if (unres_pdb) then
3182 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3183 ! print *,i,'tu dochodze'
3184 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3192 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3196 dcj=(c(j,i-2)-c(j,i-3))/2.0
3197 if (dcj.eq.0) dcj=1.23591524223
3202 else !itype(i+1,1).eq.ntyp1
3203 ! if (unres_pdb) then
3204 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3205 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3212 c(j,i)=c(j,i+1)-1.9d0*e2(j)
3216 dcj=(c(j,i+3)-c(j,i+2))/2.0
3217 if (dcj.eq.0) dcj=1.23591524223
3222 endif !itype(i+1,1).eq.ntyp1
3223 endif !itype.eq.ntyp1
3227 ! Calculate the CM of the last side chain.
3231 dc(j,ires)=sccor(j,iii)
3234 call sccenter(ires,iii,sccor)
3240 ! print *,"molecule",molecule
3241 if ((itype(nres,1).ne.10)) then
3243 if (molecule.eq.5) molecule=molecprev
3244 itype(nres,molecule)=ntyp1_molec(molecule)
3245 nres_molec(molecule)=nres_molec(molecule)+1
3246 ! if (unres_pdb) then
3247 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3248 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3255 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3259 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3260 c(j,nres)=c(j,nres-1)+dcj
3261 c(j,2*nres)=c(j,nres)
3265 ! print *,'I have',nres, nres_molec(:)
3267 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3269 if (nres.ne.nres0) then
3270 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3272 stop "Error nres value in WHAM input"
3275 !---------------------------------
3276 !el reallocate tables
3279 ! hfrag_alloc(j,i)=hfrag(j,i)
3282 ! bfrag_alloc(j,i)=bfrag(j,i)
3288 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3289 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3290 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3291 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3295 ! hfrag(j,i)=hfrag_alloc(j,i)
3300 ! bfrag(j,i)=bfrag_alloc(j,i)
3303 !el end reallocate tables
3304 !---------------------------------
3312 c(j,2*nres)=c(j,nres)
3315 if (itype(1,1).eq.ntyp1) then
3319 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3320 call refsys(2,3,4,e1,e2,e3,fail)
3327 c(j,1)=c(j,2)-1.9d0*e2(j)
3331 dcj=(c(j,4)-c(j,3))/2.0
3337 ! First lets assign correct dummy to correct type of chain
3339 if (itype(1,1).eq.ntyp1) then
3340 if (itype(2,1).eq.0) then
3342 if (itype(2,j).ne.0) then
3344 itype(1,j)=ntyp1_molec(j)
3345 nres_molec(1)=nres_molec(1)-1
3346 nres_molec(j)=nres_molec(j)+1
3353 print *,'I have',nres, nres_molec(:)
3355 ! Copy the coordinates to reference coordinates
3361 ! Calculate internal coordinates.
3363 write (iout,'(/a)') &
3364 "Cartesian coordinates of the reference structure"
3365 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3366 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3368 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3369 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3370 (c(j,ires+nres),j=1,3)
3373 ! znamy już nres wiec mozna alokowac tablice
3374 ! Calculate internal coordinates.
3375 if(me.eq.king.or..not.out1file)then
3376 write (iout,'(a)') &
3377 "Backbone and SC coordinates as read from the PDB"
3379 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3380 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3381 (c(j,nres+ires),j=1,3)
3384 ! NOW LETS ROCK! SORTING
3385 allocate(c_temporary(3,2*nres))
3386 allocate(itype_temporary(nres,5))
3387 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3388 if (.not.allocated(istype)) write(iout,*) &
3389 "SOMETHING WRONG WITH ISTYTPE"
3390 allocate(istype_temp(nres))
3391 itype_temporary(:,:)=0
3395 if (itype(i,k).ne.0) then
3397 c_temporary(j,seqalingbegin)=c(j,i)
3398 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3401 itype_temporary(seqalingbegin,k)=itype(i,k)
3402 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3403 istype_temp(seqalingbegin)=istype(i)
3404 molnum(seqalingbegin)=k
3405 seqalingbegin=seqalingbegin+1
3411 c(j,i)=c_temporary(j,i)
3416 itype(i,k)=itype_temporary(i,k)
3417 istype(i)=istype_temp(i)
3420 ! if (itype(1,1).eq.ntyp1) then
3423 ! if (unres_pdb) then
3424 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3425 ! call refsys(2,3,4,e1,e2,e3,fail)
3432 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3436 ! dcj=(c(j,4)-c(j,3))/2.0
3438 ! c(j,nres+1)=c(j,1)
3444 write (iout,'(/a)') &
3445 "Cartesian coordinates of the reference structure after sorting"
3446 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3447 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3449 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3450 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3451 (c(j,ires+nres),j=1,3)
3455 ! print *,seqalingbegin,nres
3456 if(.not.allocated(vbld)) then
3457 allocate(vbld(2*nres))
3462 if(.not.allocated(vbld_inv)) then
3463 allocate(vbld_inv(2*nres))
3469 if(.not.allocated(theta)) then
3470 allocate(theta(nres+2))
3474 if(.not.allocated(phi)) allocate(phi(nres+2))
3475 if(.not.allocated(alph)) allocate(alph(nres+2))
3476 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3477 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3478 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3479 if(.not.allocated(costtab)) allocate(costtab(nres))
3480 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3481 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3482 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3483 if(.not.allocated(xxref)) allocate(xxref(nres))
3484 if(.not.allocated(yyref)) allocate(yyref(nres))
3485 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3486 if(.not.allocated(dc_norm)) then
3487 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3488 allocate(dc_norm(3,0:2*nres+2))
3492 call int_from_cart(.true.,.false.)
3493 call sc_loc_geom(.false.)
3495 thetaref(i)=theta(i)
3505 dc(j,i)=c(j,i+1)-c(j,i)
3506 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3511 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3512 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3514 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3518 ! Copy the coordinates to reference coordinates
3519 ! Splits to single chain if occurs
3523 ! cref(j,i,cou)=c(j,i)
3527 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3528 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3529 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3530 !-----------------------------
3534 write (iout,*) "symetr", symetr
3537 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3539 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3542 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3547 cref(j,i,cou)=c(j,i)
3548 cref(j,i+nres,cou)=c(j,i+nres)
3550 chain_rep(j,lll,kkk)=c(j,i)
3551 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3555 write (iout,*) chain_length
3556 if (chain_length.eq.0) chain_length=nres
3558 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3559 chain_rep(j,chain_length+nres,symetr) &
3560 =chain_rep(j,chain_length+nres,1)
3563 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3565 ! do kkk=1,chain_length
3566 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3570 ! makes copy of chains
3571 write (iout,*) "symetr", symetr
3576 if (symetr.gt.1) then
3583 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3589 ! write (iout,*) i,icha
3590 do lll=1,chain_length
3592 if (cou.le.nres) then
3594 kupa=mod(lll,chain_length)
3595 iprzes=(kkk-1)*chain_length+lll
3596 if (kupa.eq.0) kupa=chain_length
3597 ! write (iout,*) "kupa", kupa
3598 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3599 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3606 !-koniec robienia kopii
3609 write (iout,*) "nowa struktura", nperm
3611 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3613 cref(3,i,kkk),cref(1,nres+i,kkk),&
3614 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3616 100 format (//' alpha-carbon coordinates ',&
3617 ' centroid coordinates'/ &
3618 ' ', 6X,'X',11X,'Y',11X,'Z', &
3619 10X,'X',11X,'Y',11X,'Z')
3620 110 format (a,'(',i3,')',6f12.5)
3626 bfrag(i,j)=bfrag(i,j)-ishift
3632 hfrag(i,j)=hfrag(i,j)-ishift
3638 end subroutine readpdb
3639 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3640 !-----------------------------------------------------------------------------
3642 !-----------------------------------------------------------------------------
3643 subroutine read_control
3657 use random, only: random_init
3658 ! implicit real*8 (a-h,o-z)
3659 ! include 'DIMENSIONS'
3661 use prng, only:prng_restart
3663 logical :: OKRandom!, prng_restart
3666 ! include 'COMMON.IOUNITS'
3667 ! include 'COMMON.TIME1'
3668 ! include 'COMMON.THREAD'
3669 ! include 'COMMON.SBRIDGE'
3670 ! include 'COMMON.CONTROL'
3671 ! include 'COMMON.MCM'
3672 ! include 'COMMON.MAP'
3673 ! include 'COMMON.HEADER'
3674 ! include 'COMMON.CSA'
3675 ! include 'COMMON.CHAIN'
3676 ! include 'COMMON.MUCA'
3677 ! include 'COMMON.MD'
3678 ! include 'COMMON.FFIELD'
3679 ! include 'COMMON.INTERACT'
3680 ! include 'COMMON.SETUP'
3681 !el integer :: KDIAG,ICORFL,IXDR
3682 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3683 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3684 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3685 ! character(len=80) :: ucase
3686 character(len=640) :: controlcard
3688 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3694 read (INP,'(a)') titel
3695 call card_concat(controlcard,.true.)
3696 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3697 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3698 call reada(controlcard,'SEED',seed,0.0D0)
3699 call random_init(seed)
3700 ! Set up the time limit (caution! The time must be input in minutes!)
3701 read_cart=index(controlcard,'READ_CART').gt.0
3702 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3703 call readi(controlcard,'SYM',symetr,1)
3704 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3705 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3706 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3707 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3708 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3709 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3710 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3711 call reada(controlcard,'DRMS',drms,0.1D0)
3712 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3713 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3714 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3715 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3716 write (iout,'(a,f10.1)')'DRMS = ',drms
3717 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3718 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3720 call readi(controlcard,'NZ_START',nz_start,0)
3721 call readi(controlcard,'NZ_END',nz_end,0)
3722 call readi(controlcard,'IZ_SC',iz_sc,0)
3723 timlim=60.0D0*timlim
3724 safety = 60.0d0*safety
3727 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3728 !C SHIELD keyword sets if the shielding effect of side-chains is used
3729 !C 0 denots no shielding is used all peptide are equally despite the
3730 !C solvent accesible area
3731 !C 1 the newly introduced function
3732 !C 2 reseved for further possible developement
3733 call readi(controlcard,'SHIELD',shield_mode,0)
3734 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3735 write(iout,*) "shield_mode",shield_mode
3736 !C Varibles set size of box
3737 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3738 protein=index(controlcard,"PROTEIN").gt.0
3739 ions=index(controlcard,"IONS").gt.0
3740 nucleic=index(controlcard,"NUCLEIC").gt.0
3741 write (iout,*) "with_theta_constr ",with_theta_constr
3742 AFMlog=(index(controlcard,'AFM'))
3743 selfguide=(index(controlcard,'SELFGUIDE'))
3744 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3745 call readi(controlcard,'GENCONSTR',genconstr,0)
3746 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3747 call reada(controlcard,'BOXY',boxysize,100.0d0)
3748 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3749 call readi(controlcard,'TUBEMOD',tubemode,0)
3750 print *,"SCELE",scelemode
3751 call readi(controlcard,"SCELEMODE",scelemode,0)
3752 print *,"SCELE",scelemode
3754 ! elemode = 0 is orignal UNRES electrostatics
3755 ! elemode = 1 is "Momo" potentials in progress
3756 ! elemode = 2 is in development EVALD
3757 write (iout,*) TUBEmode,"TUBEMODE"
3758 if (TUBEmode.gt.0) then
3759 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3760 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3761 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3762 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3763 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3764 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3765 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3766 buftubebot=bordtubebot+tubebufthick
3767 buftubetop=bordtubetop-tubebufthick
3770 ! CUTOFFF ON ELECTROSTATICS
3771 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3772 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3773 write(iout,*) "R_CUT_ELE=",r_cut_ele
3774 ! Lipidic parameters
3775 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3776 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3777 if (lipthick.gt.0.0d0) then
3778 bordliptop=(boxzsize+lipthick)/2.0
3779 bordlipbot=bordliptop-lipthick
3780 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3781 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3782 buflipbot=bordlipbot+lipbufthick
3783 bufliptop=bordliptop-lipbufthick
3784 if ((lipbufthick*2.0d0).gt.lipthick) &
3785 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3786 endif !lipthick.gt.0
3787 write(iout,*) "bordliptop=",bordliptop
3788 write(iout,*) "bordlipbot=",bordlipbot
3789 write(iout,*) "bufliptop=",bufliptop
3790 write(iout,*) "buflipbot=",buflipbot
3791 write (iout,*) "SHIELD MODE",shield_mode
3793 !C-------------------------
3794 minim=(index(controlcard,'MINIMIZE').gt.0)
3795 dccart=(index(controlcard,'CART').gt.0)
3796 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3797 overlapsc=.not.overlapsc
3798 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3799 searchsc=.not.searchsc
3800 sideadd=(index(controlcard,'SIDEADD').gt.0)
3801 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3802 outpdb=(index(controlcard,'PDBOUT').gt.0)
3803 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3804 pdbref=(index(controlcard,'PDBREF').gt.0)
3805 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3806 indpdb=index(controlcard,'PDBSTART')
3807 extconf=(index(controlcard,'EXTCONF').gt.0)
3808 call readi(controlcard,'IPRINT',iprint,0)
3809 call readi(controlcard,'MAXGEN',maxgen,10000)
3810 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3811 call readi(controlcard,"KDIAG",kdiag,0)
3812 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3813 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3814 write (iout,*) "RESCALE_MODE",rescale_mode
3815 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3816 if (index(controlcard,'REGULAR').gt.0.0D0) then
3817 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3821 if (index(controlcard,'CHECKGRAD').gt.0) then
3823 if (index(controlcard,'CART').gt.0) then
3825 elseif (index(controlcard,'CARINT').gt.0) then
3830 elseif (index(controlcard,'THREAD').gt.0) then
3832 call readi(controlcard,'THREAD',nthread,0)
3833 if (nthread.gt.0) then
3834 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3837 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3838 stop 'Error termination in Read_Control.'
3840 else if (index(controlcard,'MCMA').gt.0) then
3842 else if (index(controlcard,'MCEE').gt.0) then
3844 else if (index(controlcard,'MULTCONF').gt.0) then
3846 else if (index(controlcard,'MAP').gt.0) then
3848 call readi(controlcard,'MAP',nmap,0)
3849 else if (index(controlcard,'CSA').gt.0) then
3851 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3853 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3856 !fcm else if (index(controlcard,'MCMF').gt.0) then
3858 else if (index(controlcard,'SOFTREG').gt.0) then
3860 else if (index(controlcard,'CHECK_BOND').gt.0) then
3862 else if (index(controlcard,'TEST').gt.0) then
3864 else if (index(controlcard,'MD').gt.0) then
3866 else if (index(controlcard,'RE ').gt.0) then
3870 lmuca=index(controlcard,'MUCA').gt.0
3871 call readi(controlcard,'MUCADYN',mucadyn,0)
3872 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3873 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3875 write (iout,*) 'MUCADYN=',mucadyn
3876 write (iout,*) 'MUCASMOOTH=',muca_smooth
3879 iscode=index(controlcard,'ONE_LETTER')
3880 indphi=index(controlcard,'PHI')
3881 indback=index(controlcard,'BACK')
3882 iranconf=index(controlcard,'RAND_CONF')
3883 i2ndstr=index(controlcard,'USE_SEC_PRED')
3884 gradout=index(controlcard,'GRADOUT').gt.0
3885 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3886 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3887 if (me.eq.king .or. .not.out1file ) &
3888 write (iout,*) "DISTCHAINMAX",distchainmax
3890 if(me.eq.king.or..not.out1file) &
3891 write (iout,'(2a)') diagmeth(kdiag),&
3892 ' routine used to diagonalize matrices.'
3893 if (shield_mode.gt.0) then
3895 !C VSolvSphere the volume of solving sphere
3897 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3898 !C there will be no distinction between proline peptide group and normal peptide
3899 !C group in case of shielding parameters
3900 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3901 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3902 write (iout,*) VSolvSphere,VSolvSphere_div
3903 !C long axis of side chain
3905 long_r_sidechain(i)=vbldsc0(1,i)
3906 short_r_sidechain(i)=sigma0(i)
3907 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3913 end subroutine read_control
3914 !-----------------------------------------------------------------------------
3915 subroutine read_REMDpar
3917 ! Read REMD settings
3924 use control_data, only:out1file
3925 ! implicit real*8 (a-h,o-z)
3926 ! include 'DIMENSIONS'
3927 ! include 'COMMON.IOUNITS'
3928 ! include 'COMMON.TIME1'
3929 ! include 'COMMON.MD'
3932 !el include 'COMMON.LANGEVIN'
3934 !el include 'COMMON.LANGEVIN.lang0'
3936 ! include 'COMMON.INTERACT'
3937 ! include 'COMMON.NAMES'
3938 ! include 'COMMON.GEO'
3939 ! include 'COMMON.REMD'
3940 ! include 'COMMON.CONTROL'
3941 ! include 'COMMON.SETUP'
3942 ! character(len=80) :: ucase
3943 character(len=320) :: controlcard
3944 character(len=3200) :: controlcard1
3945 integer :: iremd_m_total
3948 ! real(kind=8) :: var,ene
3950 if(me.eq.king.or..not.out1file) &
3951 write (iout,*) "REMD setup"
3953 call card_concat(controlcard,.true.)
3954 call readi(controlcard,"NREP",nrep,3)
3955 call readi(controlcard,"NSTEX",nstex,1000)
3956 call reada(controlcard,"RETMIN",retmin,10.0d0)
3957 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3958 mremdsync=(index(controlcard,'SYNC').gt.0)
3959 call readi(controlcard,"NSYN",i_sync_step,100)
3960 restart1file=(index(controlcard,'REST1FILE').gt.0)
3961 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3962 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3963 if(max_cache_traj_use.gt.max_cache_traj) &
3964 max_cache_traj_use=max_cache_traj
3965 if(me.eq.king.or..not.out1file) then
3966 !d if (traj1file) then
3967 !rc caching is in testing - NTWX is not ignored
3968 !d write (iout,*) "NTWX value is ignored"
3969 !d write (iout,*) " trajectory is stored to one file by master"
3970 !d write (iout,*) " before exchange at NSTEX intervals"
3972 write (iout,*) "NREP= ",nrep
3973 write (iout,*) "NSTEX= ",nstex
3974 write (iout,*) "SYNC= ",mremdsync
3975 write (iout,*) "NSYN= ",i_sync_step
3976 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3979 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3980 if (index(controlcard,'TLIST').gt.0) then
3982 call card_concat(controlcard1,.true.)
3983 read(controlcard1,*) (remd_t(i),i=1,nrep)
3984 if(me.eq.king.or..not.out1file) &
3985 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3988 if (index(controlcard,'MLIST').gt.0) then
3990 call card_concat(controlcard1,.true.)
3991 read(controlcard1,*) (remd_m(i),i=1,nrep)
3992 if(me.eq.king.or..not.out1file) then
3993 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3996 iremd_m_total=iremd_m_total+remd_m(i)
3998 write (iout,*) 'Total number of replicas ',iremd_m_total
4001 if(me.eq.king.or..not.out1file) &
4002 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4004 end subroutine read_REMDpar
4005 !-----------------------------------------------------------------------------
4006 subroutine read_MDpar
4010 use control_data, only: r_cut,rlamb,out1file
4012 use geometry_data, only: pi
4014 ! implicit real*8 (a-h,o-z)
4015 ! include 'DIMENSIONS'
4016 ! include 'COMMON.IOUNITS'
4017 ! include 'COMMON.TIME1'
4018 ! include 'COMMON.MD'
4021 !el include 'COMMON.LANGEVIN'
4023 !el include 'COMMON.LANGEVIN.lang0'
4025 ! include 'COMMON.INTERACT'
4026 ! include 'COMMON.NAMES'
4027 ! include 'COMMON.GEO'
4028 ! include 'COMMON.SETUP'
4029 ! include 'COMMON.CONTROL'
4030 ! include 'COMMON.SPLITELE'
4031 ! character(len=80) :: ucase
4032 character(len=320) :: controlcard
4037 call card_concat(controlcard,.true.)
4038 call readi(controlcard,"NSTEP",n_timestep,1000000)
4039 call readi(controlcard,"NTWE",ntwe,100)
4040 call readi(controlcard,"NTWX",ntwx,1000)
4041 call reada(controlcard,"DT",d_time,1.0d-1)
4042 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4043 call reada(controlcard,"DAMAX",damax,1.0d1)
4044 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4045 call readi(controlcard,"LANG",lang,0)
4046 RESPA = index(controlcard,"RESPA") .gt. 0
4047 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4048 ntime_split0=ntime_split
4049 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4050 ntime_split0=ntime_split
4051 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4052 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4053 rest = index(controlcard,"REST").gt.0
4054 tbf = index(controlcard,"TBF").gt.0
4055 usampl = index(controlcard,"USAMPL").gt.0
4056 mdpdb = index(controlcard,"MDPDB").gt.0
4057 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4058 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4059 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4060 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4061 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4062 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4063 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4064 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4065 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4066 large = index(controlcard,"LARGE").gt.0
4067 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4068 rattle = index(controlcard,"RATTLE").gt.0
4069 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4075 if(me.eq.king.or..not.out1file) then
4077 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4079 write (iout,'(a)') "The units are:"
4080 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4081 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4082 " acceleration: angstrom/(48.9 fs)**2"
4083 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4085 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4086 write (iout,'(a60,f10.5,a)') &
4087 "Initial time step of numerical integration:",d_time,&
4089 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4091 write (iout,'(2a,i4,a)') &
4092 "A-MTS algorithm used; initial time step for fast-varying",&
4093 " short-range forces split into",ntime_split," steps."
4094 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4095 r_cut," lambda",rlamb
4097 write (iout,'(2a,f10.5)') &
4098 "Maximum acceleration threshold to reduce the time step",&
4099 "/increase split number:",damax
4100 write (iout,'(2a,f10.5)') &
4101 "Maximum predicted energy drift to reduce the timestep",&
4102 "/increase split number:",edriftmax
4103 write (iout,'(a60,f10.5)') &
4104 "Maximum velocity threshold to reduce velocities:",dvmax
4105 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4106 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4107 if (rattle) write (iout,'(a60)') &
4108 "Rattle algorithm used to constrain the virtual bonds"
4112 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4113 call reada(controlcard,"RWAT",rwat,1.4d0)
4114 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4115 surfarea=index(controlcard,"SURFAREA").gt.0
4116 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4117 if(me.eq.king.or..not.out1file)then
4118 write (iout,'(/a,$)') "Langevin dynamics calculation"
4120 write (iout,'(a/)') &
4121 " with direct integration of Langevin equations"
4122 else if (lang.eq.2) then
4123 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4124 else if (lang.eq.3) then
4125 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4126 else if (lang.eq.4) then
4127 write (iout,'(a/)') " in overdamped mode"
4129 write (iout,'(//a,i5)') &
4130 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4133 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4134 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4135 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4136 write (iout,'(a60,f10.5)') &
4137 "Scaling factor of the friction forces:",scal_fric
4138 if (surfarea) write (iout,'(2a,i10,a)') &
4139 "Friction coefficients will be scaled by solvent-accessible",&
4140 " surface area every",reset_fricmat," steps."
4142 ! Calculate friction coefficients and bounds of stochastic forces
4143 eta=6*pi*cPoise*etawat
4144 if(me.eq.king.or..not.out1file) &
4145 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4148 do j=1,5 !types of molecules
4149 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4150 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4152 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4153 do j=1,5 !types of molecules
4155 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4156 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4160 if(me.eq.king.or..not.out1file)then
4161 write (iout,'(/2a/)') &
4162 "Radii of site types and friction coefficients and std's of",&
4163 " stochastic forces of fully exposed sites"
4164 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4166 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4167 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4171 if(me.eq.king.or..not.out1file)then
4172 write (iout,'(a)') "Berendsen bath calculation"
4173 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4174 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4176 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4177 count_reset_moment," steps"
4179 write (iout,'(a,i10,a)') &
4180 "Velocities will be reset at random every",count_reset_vel,&
4184 if(me.eq.king.or..not.out1file) &
4185 write (iout,'(a31)') "Microcanonical mode calculation"
4187 if(me.eq.king.or..not.out1file)then
4188 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4190 write(iout,*) "MD running with constraints."
4191 write(iout,*) "Equilibration time ", eq_time, " mtus."
4192 write(iout,*) "Constraining ", nfrag," fragments."
4193 write(iout,*) "Length of each fragment, weight and q0:"
4195 write (iout,*) "Set of restraints #",iset
4197 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4198 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4200 write(iout,*) "constraints between ", npair, "fragments."
4201 write(iout,*) "constraint pairs, weights and q0:"
4203 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4204 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4206 write(iout,*) "angle constraints within ", nfrag_back,&
4207 "backbone fragments."
4208 write(iout,*) "fragment, weights:"
4210 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4211 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4212 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4215 iset=mod(kolor,nset)+1
4218 if(me.eq.king.or..not.out1file) &
4219 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4221 end subroutine read_MDpar
4222 !-----------------------------------------------------------------------------
4226 ! implicit real*8 (a-h,o-z)
4227 ! include 'DIMENSIONS'
4228 ! include 'COMMON.MAP'
4229 ! include 'COMMON.IOUNITS'
4230 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4231 character(len=80) :: mapcard !,ucase
4234 ! real(kind=8) :: var,ene
4237 read (inp,'(a)') mapcard
4238 mapcard=ucase(mapcard)
4239 if (index(mapcard,'PHI').gt.0) then
4241 else if (index(mapcard,'THE').gt.0) then
4243 else if (index(mapcard,'ALP').gt.0) then
4245 else if (index(mapcard,'OME').gt.0) then
4248 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4249 stop 'Error - illegal variable spec in MAP card.'
4251 call readi (mapcard,'RES1',res1(imap),0)
4252 call readi (mapcard,'RES2',res2(imap),0)
4253 if (res1(imap).eq.0) then
4254 res1(imap)=res2(imap)
4255 else if (res2(imap).eq.0) then
4256 res2(imap)=res1(imap)
4258 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4259 write (iout,'(a)') &
4260 'Error - illegal definition of variable group in MAP.'
4261 stop 'Error - illegal definition of variable group in MAP.'
4263 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4264 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4265 call readi(mapcard,'NSTEP',nstep(imap),0)
4266 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4267 write (iout,'(a)') &
4268 'Illegal boundary and/or step size specification in MAP.'
4269 stop 'Illegal boundary and/or step size specification in MAP.'
4273 end subroutine map_read
4274 !-----------------------------------------------------------------------------
4277 use control_data, only: vdisulf
4279 ! implicit real*8 (a-h,o-z)
4280 ! include 'DIMENSIONS'
4281 ! include 'COMMON.IOUNITS'
4282 ! include 'COMMON.GEO'
4283 ! include 'COMMON.CSA'
4284 ! include 'COMMON.BANK'
4285 ! include 'COMMON.CONTROL'
4286 ! character(len=80) :: ucase
4287 character(len=620) :: mcmcard
4289 ! integer :: ntf,ik,iw_pdb
4290 ! real(kind=8) :: var,ene
4292 call card_concat(mcmcard,.true.)
4294 call readi(mcmcard,'NCONF',nconf,50)
4295 call readi(mcmcard,'NADD',nadd,0)
4296 call readi(mcmcard,'JSTART',jstart,1)
4297 call readi(mcmcard,'JEND',jend,1)
4298 call readi(mcmcard,'NSTMAX',nstmax,500000)
4299 call readi(mcmcard,'N0',n0,1)
4300 call readi(mcmcard,'N1',n1,6)
4301 call readi(mcmcard,'N2',n2,4)
4302 call readi(mcmcard,'N3',n3,0)
4303 call readi(mcmcard,'N4',n4,0)
4304 call readi(mcmcard,'N5',n5,0)
4305 call readi(mcmcard,'N6',n6,10)
4306 call readi(mcmcard,'N7',n7,0)
4307 call readi(mcmcard,'N8',n8,0)
4308 call readi(mcmcard,'N9',n9,0)
4309 call readi(mcmcard,'N14',n14,0)
4310 call readi(mcmcard,'N15',n15,0)
4311 call readi(mcmcard,'N16',n16,0)
4312 call readi(mcmcard,'N17',n17,0)
4313 call readi(mcmcard,'N18',n18,0)
4315 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4317 call readi(mcmcard,'NDIFF',ndiff,2)
4318 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4319 call readi(mcmcard,'IS1',is1,1)
4320 call readi(mcmcard,'IS2',is2,8)
4321 call readi(mcmcard,'NRAN0',nran0,4)
4322 call readi(mcmcard,'NRAN1',nran1,2)
4323 call readi(mcmcard,'IRR',irr,1)
4324 call readi(mcmcard,'NSEED',nseed,20)
4325 call readi(mcmcard,'NTOTAL',ntotal,10000)
4326 call reada(mcmcard,'CUT1',cut1,2.0d0)
4327 call reada(mcmcard,'CUT2',cut2,5.0d0)
4328 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4329 call readi(mcmcard,'ICMAX',icmax,3)
4330 call readi(mcmcard,'IRESTART',irestart,0)
4331 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4334 call reada(mcmcard,'DELE',dele,20.0d0)
4335 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4336 call readi(mcmcard,'IREF',iref,0)
4337 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4338 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4339 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4340 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4341 write (iout,*) "NCONF_IN",nconf_in
4343 end subroutine csaread
4344 !-----------------------------------------------------------------------------
4348 use control_data, only: MaxMoveType
4351 ! implicit real*8 (a-h,o-z)
4352 ! include 'DIMENSIONS'
4353 ! include 'COMMON.MCM'
4354 ! include 'COMMON.MCE'
4355 ! include 'COMMON.IOUNITS'
4356 ! character(len=80) :: ucase
4357 character(len=320) :: mcmcard
4360 ! real(kind=8) :: var,ene
4362 call card_concat(mcmcard,.true.)
4363 call readi(mcmcard,'MAXACC',maxacc,100)
4364 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4365 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4366 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4367 call readi(mcmcard,'MAXREPM',maxrepm,200)
4368 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4369 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4370 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4371 call reada(mcmcard,'E_UP',e_up,5.0D0)
4372 call reada(mcmcard,'DELTE',delte,0.1D0)
4373 call readi(mcmcard,'NSWEEP',nsweep,5)
4374 call readi(mcmcard,'NSTEPH',nsteph,0)
4375 call readi(mcmcard,'NSTEPC',nstepc,0)
4376 call reada(mcmcard,'TMIN',tmin,298.0D0)
4377 call reada(mcmcard,'TMAX',tmax,298.0D0)
4378 call readi(mcmcard,'NWINDOW',nwindow,0)
4379 call readi(mcmcard,'PRINT_MC',print_mc,0)
4380 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4381 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4382 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4383 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4384 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4385 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4386 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4387 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4388 if (nwindow.gt.0) then
4389 allocate(winstart(nwindow)) !!el (maxres)
4390 allocate(winend(nwindow)) !!el
4391 allocate(winlen(nwindow)) !!el
4392 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4394 winlen(i)=winend(i)-winstart(i)+1
4397 if (tmax.lt.tmin) tmax=tmin
4398 if (tmax.eq.tmin) then
4402 if (nstepc.gt.0 .and. nsteph.gt.0) then
4403 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4404 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4406 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4407 ! Probabilities of different move types
4408 sumpro_type(0)=0.0D0
4409 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4410 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4411 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4412 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4413 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4414 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4415 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4417 print *,'i',i,' sumprotype',sumpro_type(i)
4418 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4419 print *,'i',i,' sumprotype',sumpro_type(i)
4422 end subroutine mcmread
4423 !-----------------------------------------------------------------------------
4424 subroutine read_minim
4427 ! implicit real*8 (a-h,o-z)
4428 ! include 'DIMENSIONS'
4429 ! include 'COMMON.MINIM'
4430 ! include 'COMMON.IOUNITS'
4431 ! character(len=80) :: ucase
4432 character(len=320) :: minimcard
4434 ! integer :: ntf,ik,iw_pdb
4435 ! real(kind=8) :: var,ene
4437 call card_concat(minimcard,.true.)
4438 call readi(minimcard,'MAXMIN',maxmin,2000)
4439 call readi(minimcard,'MAXFUN',maxfun,5000)
4440 call readi(minimcard,'MINMIN',minmin,maxmin)
4441 call readi(minimcard,'MINFUN',minfun,maxmin)
4442 call reada(minimcard,'TOLF',tolf,1.0D-2)
4443 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4444 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4445 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4446 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4447 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4448 'Options in energy minimization:'
4449 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4450 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4451 'MinMin:',MinMin,' MinFun:',MinFun,&
4452 ' TolF:',TolF,' RTolF:',RTolF
4454 end subroutine read_minim
4455 !-----------------------------------------------------------------------------
4456 subroutine openunits
4458 use MD_data, only: usampl
4461 use control_data, only:out1file
4462 use control, only: getenv_loc
4463 ! implicit real*8 (a-h,o-z)
4464 ! include 'DIMENSIONS'
4467 character(len=16) :: form,nodename
4468 integer :: nodelen,ierror,npos
4470 ! include 'COMMON.SETUP'
4471 ! include 'COMMON.IOUNITS'
4472 ! include 'COMMON.MD'
4473 ! include 'COMMON.CONTROL'
4474 integer :: lenpre,lenpot,lentmp !,ilen
4476 character(len=3) :: out1file_text !,ucase
4477 character(len=3) :: ll
4480 ! integer :: ntf,ik,iw_pdb
4481 ! real(kind=8) :: var,ene
4483 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4484 call getenv_loc("PREFIX",prefix)
4486 call getenv_loc("POT",pot)
4487 call getenv_loc("DIRTMP",tmpdir)
4488 call getenv_loc("CURDIR",curdir)
4489 call getenv_loc("OUT1FILE",out1file_text)
4490 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4491 out1file_text=ucase(out1file_text)
4492 if (out1file_text(1:1).eq."Y") then
4495 out1file=fg_rank.gt.0
4500 if (lentmp.gt.0) then
4501 write (*,'(80(1h!))')
4502 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4503 write (*,'(80(1h!))')
4504 write (*,*)"All output files will be on node /tmp directory."
4506 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4507 if (me.eq.king) then
4508 write (*,*) "The master node is ",nodename
4509 else if (fg_rank.eq.0) then
4510 write (*,*) "I am the CG slave node ",nodename
4512 write (*,*) "I am the FG slave node ",nodename
4515 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4516 lenpre = lentmp+lenpre+1
4518 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4519 ! Get the names and open the input files
4520 #if defined(WINIFL) || defined(WINPGI)
4521 open(1,file=pref_orig(:ilen(pref_orig))// &
4522 '.inp',status='old',readonly,shared)
4523 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4524 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4525 ! Get parameter filenames and open the parameter files.
4526 call getenv_loc('BONDPAR',bondname)
4527 open (ibond,file=bondname,status='old',readonly,shared)
4528 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4529 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4530 call getenv_loc('THETPAR',thetname)
4531 open (ithep,file=thetname,status='old',readonly,shared)
4532 call getenv_loc('ROTPAR',rotname)
4533 open (irotam,file=rotname,status='old',readonly,shared)
4534 call getenv_loc('TORPAR',torname)
4535 open (itorp,file=torname,status='old',readonly,shared)
4536 call getenv_loc('TORDPAR',tordname)
4537 open (itordp,file=tordname,status='old',readonly,shared)
4538 call getenv_loc('FOURIER',fouriername)
4539 open (ifourier,file=fouriername,status='old',readonly,shared)
4540 call getenv_loc('ELEPAR',elename)
4541 open (ielep,file=elename,status='old',readonly,shared)
4542 call getenv_loc('SIDEPAR',sidename)
4543 open (isidep,file=sidename,status='old',readonly,shared)
4545 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4546 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4547 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4548 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4549 call getenv_loc('TORPAR_NUCL',torname_nucl)
4550 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4551 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4552 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4553 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4554 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4557 #elif (defined CRAY) || (defined AIX)
4558 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4560 ! print *,"Processor",myrank," opened file 1"
4561 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4562 ! print *,"Processor",myrank," opened file 9"
4563 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4564 ! Get parameter filenames and open the parameter files.
4565 call getenv_loc('BONDPAR',bondname)
4566 open (ibond,file=bondname,status='old',action='read')
4567 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4568 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4570 ! print *,"Processor",myrank," opened file IBOND"
4571 call getenv_loc('THETPAR',thetname)
4572 open (ithep,file=thetname,status='old',action='read')
4573 ! print *,"Processor",myrank," opened file ITHEP"
4574 call getenv_loc('ROTPAR',rotname)
4575 open (irotam,file=rotname,status='old',action='read')
4576 ! print *,"Processor",myrank," opened file IROTAM"
4577 call getenv_loc('TORPAR',torname)
4578 open (itorp,file=torname,status='old',action='read')
4579 ! print *,"Processor",myrank," opened file ITORP"
4580 call getenv_loc('TORDPAR',tordname)
4581 open (itordp,file=tordname,status='old',action='read')
4582 ! print *,"Processor",myrank," opened file ITORDP"
4583 call getenv_loc('SCCORPAR',sccorname)
4584 open (isccor,file=sccorname,status='old',action='read')
4585 ! print *,"Processor",myrank," opened file ISCCOR"
4586 call getenv_loc('FOURIER',fouriername)
4587 open (ifourier,file=fouriername,status='old',action='read')
4588 ! print *,"Processor",myrank," opened file IFOURIER"
4589 call getenv_loc('ELEPAR',elename)
4590 open (ielep,file=elename,status='old',action='read')
4591 ! print *,"Processor",myrank," opened file IELEP"
4592 call getenv_loc('SIDEPAR',sidename)
4593 open (isidep,file=sidename,status='old',action='read')
4595 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4596 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4597 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4598 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4599 call getenv_loc('TORPAR_NUCL',torname_nucl)
4600 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4601 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4602 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4603 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4604 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4606 call getenv_loc('LIPTRANPAR',liptranname)
4607 open (iliptranpar,file=liptranname,status='old',action='read')
4608 call getenv_loc('TUBEPAR',tubename)
4609 open (itube,file=tubename,status='old',action='read')
4610 call getenv_loc('IONPAR',ionname)
4611 open (iion,file=ionname,status='old',action='read')
4613 ! print *,"Processor",myrank," opened file ISIDEP"
4614 ! print *,"Processor",myrank," opened parameter files"
4616 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4617 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4618 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4619 ! Get parameter filenames and open the parameter files.
4620 call getenv_loc('BONDPAR',bondname)
4621 open (ibond,file=bondname,status='old')
4622 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4623 open (ibond_nucl,file=bondname_nucl,status='old')
4625 call getenv_loc('THETPAR',thetname)
4626 open (ithep,file=thetname,status='old')
4627 call getenv_loc('ROTPAR',rotname)
4628 open (irotam,file=rotname,status='old')
4629 call getenv_loc('TORPAR',torname)
4630 open (itorp,file=torname,status='old')
4631 call getenv_loc('TORDPAR',tordname)
4632 open (itordp,file=tordname,status='old')
4633 call getenv_loc('SCCORPAR',sccorname)
4634 open (isccor,file=sccorname,status='old')
4635 call getenv_loc('FOURIER',fouriername)
4636 open (ifourier,file=fouriername,status='old')
4637 call getenv_loc('ELEPAR',elename)
4638 open (ielep,file=elename,status='old')
4639 call getenv_loc('SIDEPAR',sidename)
4640 open (isidep,file=sidename,status='old')
4642 open (ithep_nucl,file=thetname_nucl,status='old')
4643 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4644 open (irotam_nucl,file=rotname_nucl,status='old')
4645 call getenv_loc('TORPAR_NUCL',torname_nucl)
4646 open (itorp_nucl,file=torname_nucl,status='old')
4647 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4648 open (itordp_nucl,file=tordname_nucl,status='old')
4649 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4650 open (isidep_nucl,file=sidename_nucl,status='old')
4652 call getenv_loc('LIPTRANPAR',liptranname)
4653 open (iliptranpar,file=liptranname,status='old')
4654 call getenv_loc('TUBEPAR',tubename)
4655 open (itube,file=tubename,status='old')
4656 call getenv_loc('IONPAR',ionname)
4657 open (iion,file=ionname,status='old')
4659 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4661 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4662 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4663 ! Get parameter filenames and open the parameter files.
4664 call getenv_loc('BONDPAR',bondname)
4665 open (ibond,file=bondname,status='old',action='read')
4666 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4667 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4668 call getenv_loc('THETPAR',thetname)
4669 open (ithep,file=thetname,status='old',action='read')
4670 call getenv_loc('ROTPAR',rotname)
4671 open (irotam,file=rotname,status='old',action='read')
4672 call getenv_loc('TORPAR',torname)
4673 open (itorp,file=torname,status='old',action='read')
4674 call getenv_loc('TORDPAR',tordname)
4675 open (itordp,file=tordname,status='old',action='read')
4676 call getenv_loc('SCCORPAR',sccorname)
4677 open (isccor,file=sccorname,status='old',action='read')
4679 call getenv_loc('THETPARPDB',thetname_pdb)
4680 print *,"thetname_pdb ",thetname_pdb
4681 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4682 print *,ithep_pdb," opened"
4684 call getenv_loc('FOURIER',fouriername)
4685 open (ifourier,file=fouriername,status='old',readonly)
4686 call getenv_loc('ELEPAR',elename)
4687 open (ielep,file=elename,status='old',readonly)
4688 call getenv_loc('SIDEPAR',sidename)
4689 open (isidep,file=sidename,status='old',readonly)
4691 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4692 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4693 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4694 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4695 call getenv_loc('TORPAR_NUCL',torname_nucl)
4696 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4697 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4698 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4699 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4700 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4701 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4702 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4703 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4704 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4705 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4706 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4707 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4708 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4711 call getenv_loc('LIPTRANPAR',liptranname)
4712 open (iliptranpar,file=liptranname,status='old',action='read')
4713 call getenv_loc('TUBEPAR',tubename)
4714 open (itube,file=tubename,status='old',action='read')
4715 call getenv_loc('IONPAR',ionname)
4716 open (iion,file=ionname,status='old',action='read')
4719 call getenv_loc('ROTPARPDB',rotname_pdb)
4720 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4723 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4724 #if defined(WINIFL) || defined(WINPGI)
4725 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4726 #elif (defined CRAY) || (defined AIX)
4727 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4729 open (iscpp_nucl,file=scpname_nucl,status='old')
4731 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4736 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4737 ! Use -DOLDSCP to use hard-coded constants instead.
4739 call getenv_loc('SCPPAR',scpname)
4740 #if defined(WINIFL) || defined(WINPGI)
4741 open (iscpp,file=scpname,status='old',readonly,shared)
4742 #elif (defined CRAY) || (defined AIX)
4743 open (iscpp,file=scpname,status='old',action='read')
4745 open (iscpp,file=scpname,status='old')
4747 open (iscpp,file=scpname,status='old',action='read')
4750 call getenv_loc('PATTERN',patname)
4751 #if defined(WINIFL) || defined(WINPGI)
4752 open (icbase,file=patname,status='old',readonly,shared)
4753 #elif (defined CRAY) || (defined AIX)
4754 open (icbase,file=patname,status='old',action='read')
4756 open (icbase,file=patname,status='old')
4758 open (icbase,file=patname,status='old',action='read')
4761 ! Open output file only for CG processes
4762 ! print *,"Processor",myrank," fg_rank",fg_rank
4763 if (fg_rank.eq.0) then
4765 if (nodes.eq.1) then
4768 npos = dlog10(dfloat(nodes-1))+1
4770 if (npos.lt.3) npos=3
4771 write (liczba,'(i1)') npos
4772 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4774 write (liczba,form) me
4775 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4776 liczba(:ilen(liczba))
4777 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4779 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4781 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4782 liczba(:ilen(liczba))//'.mol2'
4783 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4784 liczba(:ilen(liczba))//'.stat'
4786 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4787 //liczba(:ilen(liczba))//'.stat')
4788 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4791 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4792 liczba(:ilen(liczba))//'.const'
4797 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4798 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4799 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4800 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4801 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4803 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4805 rest2name=prefix(:ilen(prefix))//'.rst'
4807 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4810 #if defined(AIX) || defined(PGI)
4811 if (me.eq.king .or. .not. out1file) &
4812 open(iout,file=outname,status='unknown')
4814 if (fg_rank.gt.0) then
4815 write (liczba,'(i3.3)') myrank/nfgtasks
4816 write (ll,'(bz,i3.3)') fg_rank
4817 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4822 open(igeom,file=intname,status='unknown',position='append')
4823 open(ipdb,file=pdbname,status='unknown')
4824 open(imol2,file=mol2name,status='unknown')
4825 open(istat,file=statname,status='unknown',position='append')
4827 !1out open(iout,file=outname,status='unknown')
4830 if (me.eq.king .or. .not.out1file) &
4831 open(iout,file=outname,status='unknown')
4833 if (fg_rank.gt.0) then
4834 write (liczba,'(i3.3)') myrank/nfgtasks
4835 write (ll,'(bz,i3.3)') fg_rank
4836 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4841 open(igeom,file=intname,status='unknown',access='append')
4842 open(ipdb,file=pdbname,status='unknown')
4843 open(imol2,file=mol2name,status='unknown')
4844 open(istat,file=statname,status='unknown',access='append')
4846 !1out open(iout,file=outname,status='unknown')
4849 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4850 csa_seed=prefix(:lenpre)//'.CSA.seed'
4851 csa_history=prefix(:lenpre)//'.CSA.history'
4852 csa_bank=prefix(:lenpre)//'.CSA.bank'
4853 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4854 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4855 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4856 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4857 csa_int=prefix(:lenpre)//'.int'
4858 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4859 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4860 csa_in=prefix(:lenpre)//'.CSA.in'
4861 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4864 write (iout,'(80(1h-))')
4865 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4866 write (iout,'(80(1h-))')
4867 write (iout,*) "Input file : ",&
4868 pref_orig(:ilen(pref_orig))//'.inp'
4869 write (iout,*) "Output file : ",&
4870 outname(:ilen(outname))
4872 write (iout,*) "Sidechain potential file : ",&
4873 sidename(:ilen(sidename))
4875 write (iout,*) "SCp potential file : ",&
4876 scpname(:ilen(scpname))
4878 write (iout,*) "Electrostatic potential file : ",&
4879 elename(:ilen(elename))
4880 write (iout,*) "Cumulant coefficient file : ",&
4881 fouriername(:ilen(fouriername))
4882 write (iout,*) "Torsional parameter file : ",&
4883 torname(:ilen(torname))
4884 write (iout,*) "Double torsional parameter file : ",&
4885 tordname(:ilen(tordname))
4886 write (iout,*) "SCCOR parameter file : ",&
4887 sccorname(:ilen(sccorname))
4888 write (iout,*) "Bond & inertia constant file : ",&
4889 bondname(:ilen(bondname))
4890 write (iout,*) "Bending parameter file : ",&
4891 thetname(:ilen(thetname))
4892 write (iout,*) "Rotamer parameter file : ",&
4893 rotname(:ilen(rotname))
4896 write (iout,*) "Thetpdb parameter file : ",&
4897 thetname_pdb(:ilen(thetname_pdb))
4900 write (iout,*) "Threading database : ",&
4901 patname(:ilen(patname))
4903 write (iout,*)" DIRTMP : ",&
4905 write (iout,'(80(1h-))')
4908 end subroutine openunits
4909 !-----------------------------------------------------------------------------
4912 use geometry_data, only: nres,dc
4914 ! implicit real*8 (a-h,o-z)
4915 ! include 'DIMENSIONS'
4916 ! include 'COMMON.CHAIN'
4917 ! include 'COMMON.IOUNITS'
4918 ! include 'COMMON.MD'
4921 ! real(kind=8) :: var,ene
4923 open(irest2,file=rest2name,status='unknown')
4924 read(irest2,*) totT,EK,potE,totE,t_bath
4927 ! AL 4/17/17: Now reading d_t(0,:) too
4929 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4932 ! AL 4/17/17: Now reading d_c(0,:) too
4934 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4937 read (irest2,*) iset
4941 end subroutine readrst
4942 !-----------------------------------------------------------------------------
4943 subroutine read_fragments
4947 use control_data, only:out1file
4950 ! implicit real*8 (a-h,o-z)
4951 ! include 'DIMENSIONS'
4955 ! include 'COMMON.SETUP'
4956 ! include 'COMMON.CHAIN'
4957 ! include 'COMMON.IOUNITS'
4958 ! include 'COMMON.MD'
4959 ! include 'COMMON.CONTROL'
4962 ! real(kind=8) :: var,ene
4964 read(inp,*) nset,nfrag,npair,nfrag_back
4966 !el from module energy
4967 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4968 if(.not.allocated(wfrag_back)) then
4969 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4970 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4972 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4973 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4975 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4976 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4979 if(me.eq.king.or..not.out1file) &
4980 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4981 " nfrag_back",nfrag_back
4983 read(inp,*) mset(iset)
4985 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4987 if(me.eq.king.or..not.out1file) &
4988 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4989 ifrag(2,i,iset), qinfrag(i,iset)
4992 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4994 if(me.eq.king.or..not.out1file) &
4995 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4996 ipair(2,i,iset), qinpair(i,iset)
4999 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5000 wfrag_back(3,i,iset),&
5001 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5002 if(me.eq.king.or..not.out1file) &
5003 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5004 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5008 end subroutine read_fragments
5009 !-----------------------------------------------------------------------------
5011 !-----------------------------------------------------------------------------
5015 ! implicit real*8 (a-h,o-z)
5016 ! include 'DIMENSIONS'
5017 ! include 'COMMON.CSA'
5018 ! include 'COMMON.BANK'
5019 ! include 'COMMON.IOUNITS'
5021 ! integer :: ntf,ik,iw_pdb
5022 ! real(kind=8) :: var,ene
5024 open(icsa_in,file=csa_in,status="old",err=100)
5025 read(icsa_in,*) nconf
5026 read(icsa_in,*) jstart,jend
5027 read(icsa_in,*) nstmax
5028 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5029 read(icsa_in,*) nran0,nran1,irr
5030 read(icsa_in,*) nseed
5031 read(icsa_in,*) ntotal,cut1,cut2
5032 read(icsa_in,*) estop
5033 read(icsa_in,*) icmax,irestart
5034 read(icsa_in,*) ntbankm,dele,difcut
5035 read(icsa_in,*) iref,rmscut,pnccut
5036 read(icsa_in,*) ndiff
5043 end subroutine csa_read
5044 !-----------------------------------------------------------------------------
5045 subroutine initial_write
5048 ! implicit real*8 (a-h,o-z)
5049 ! include 'DIMENSIONS'
5050 ! include 'COMMON.CSA'
5051 ! include 'COMMON.BANK'
5052 ! include 'COMMON.IOUNITS'
5054 ! integer :: ntf,ik,iw_pdb
5055 ! real(kind=8) :: var,ene
5057 open(icsa_seed,file=csa_seed,status="unknown")
5058 write(icsa_seed,*) "seed"
5060 #if defined(AIX) || defined(PGI)
5061 open(icsa_history,file=csa_history,status="unknown",&
5064 open(icsa_history,file=csa_history,status="unknown",&
5067 write(icsa_history,*) nconf
5068 write(icsa_history,*) jstart,jend
5069 write(icsa_history,*) nstmax
5070 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5071 write(icsa_history,*) nran0,nran1,irr
5072 write(icsa_history,*) nseed
5073 write(icsa_history,*) ntotal,cut1,cut2
5074 write(icsa_history,*) estop
5075 write(icsa_history,*) icmax,irestart
5076 write(icsa_history,*) ntbankm,dele,difcut
5077 write(icsa_history,*) iref,rmscut,pnccut
5078 write(icsa_history,*) ndiff
5080 write(icsa_history,*)
5083 open(icsa_bank1,file=csa_bank1,status="unknown")
5084 write(icsa_bank1,*) 0
5088 end subroutine initial_write
5089 !-----------------------------------------------------------------------------
5090 subroutine restart_write
5093 ! implicit real*8 (a-h,o-z)
5094 ! include 'DIMENSIONS'
5095 ! include 'COMMON.IOUNITS'
5096 ! include 'COMMON.CSA'
5097 ! include 'COMMON.BANK'
5099 ! integer :: ntf,ik,iw_pdb
5100 ! real(kind=8) :: var,ene
5102 #if defined(AIX) || defined(PGI)
5103 open(icsa_history,file=csa_history,position="append")
5105 open(icsa_history,file=csa_history,access="append")
5107 write(icsa_history,*)
5108 write(icsa_history,*) "This is restart"
5109 write(icsa_history,*)
5110 write(icsa_history,*) nconf
5111 write(icsa_history,*) jstart,jend
5112 write(icsa_history,*) nstmax
5113 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5114 write(icsa_history,*) nran0,nran1,irr
5115 write(icsa_history,*) nseed
5116 write(icsa_history,*) ntotal,cut1,cut2
5117 write(icsa_history,*) estop
5118 write(icsa_history,*) icmax,irestart
5119 write(icsa_history,*) ntbankm,dele,difcut
5120 write(icsa_history,*) iref,rmscut,pnccut
5121 write(icsa_history,*) ndiff
5122 write(icsa_history,*)
5123 write(icsa_history,*) "irestart is: ", irestart
5125 write(icsa_history,*)
5129 end subroutine restart_write
5130 !-----------------------------------------------------------------------------
5132 !-----------------------------------------------------------------------------
5133 subroutine write_pdb(npdb,titelloc,ee)
5135 ! implicit real*8 (a-h,o-z)
5136 ! include 'DIMENSIONS'
5137 ! include 'COMMON.IOUNITS'
5138 character(len=50) :: titelloc1
5139 character*(*) :: titelloc
5140 character(len=3) :: zahl
5141 character(len=5) :: liczba5
5143 integer :: npdb !,ilen
5147 ! real(kind=8) :: var,ene
5151 if (npdb.lt.1000) then
5152 call numstr(npdb,zahl)
5153 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5155 if (npdb.lt.10000) then
5156 write(liczba5,'(i1,i4)') 0,npdb
5158 write(liczba5,'(i5)') npdb
5160 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5162 call pdbout(ee,titelloc1,ipdb)
5165 end subroutine write_pdb
5166 !-----------------------------------------------------------------------------
5168 !-----------------------------------------------------------------------------
5169 subroutine write_thread_summary
5170 ! Thread the sequence through a database of known structures
5171 use control_data, only: refstr
5173 use energy_data, only: n_ene_comp
5175 ! implicit real*8 (a-h,o-z)
5176 ! include 'DIMENSIONS'
5178 use MPI_data !include 'COMMON.INFO'
5181 ! include 'COMMON.CONTROL'
5182 ! include 'COMMON.CHAIN'
5183 ! include 'COMMON.DBASE'
5184 ! include 'COMMON.INTERACT'
5185 ! include 'COMMON.VAR'
5186 ! include 'COMMON.THREAD'
5187 ! include 'COMMON.FFIELD'
5188 ! include 'COMMON.SBRIDGE'
5189 ! include 'COMMON.HEADER'
5190 ! include 'COMMON.NAMES'
5191 ! include 'COMMON.IOUNITS'
5192 ! include 'COMMON.TIME1'
5194 integer,dimension(maxthread) :: ip
5195 real(kind=8),dimension(0:n_ene) :: energia
5197 integer :: i,j,ii,jj,ipj,ik,kk,ist
5198 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5200 write (iout,'(30x,a/)') &
5201 ' *********** Summary threading statistics ************'
5202 write (iout,'(a)') 'Initial energies:'
5203 write (iout,'(a4,2x,a12,14a14,3a8)') &
5204 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5205 'RMSnat','NatCONT','NNCONT','RMS'
5206 ! Energy sort patterns
5211 enet=ener(n_ene-1,ip(i))
5214 if (ener(n_ene-1,ip(j)).lt.enet) then
5216 enet=ener(n_ene-1,ip(j))
5228 ist=nres_base(2,ii)+ipatt(2,i)
5230 energia(i)=ener0(kk,i)
5232 etot=ener0(n_ene_comp+1,i)
5233 rmsnat=ener0(n_ene_comp+2,i)
5234 rms=ener0(n_ene_comp+3,i)
5235 frac=ener0(n_ene_comp+4,i)
5236 frac_nn=ener0(n_ene_comp+5,i)
5239 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5240 i,str_nam(ii),ist+1,&
5241 (energia(print_order(kk)),kk=1,nprint_ene),&
5242 etot,rmsnat,frac,frac_nn,rms
5244 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5245 i,str_nam(ii),ist+1,&
5246 (energia(print_order(kk)),kk=1,nprint_ene),etot
5249 write (iout,'(//a)') 'Final energies:'
5250 write (iout,'(a4,2x,a12,17a14,3a8)') &
5251 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5252 'RMSnat','NatCONT','NNCONT','RMS'
5256 ist=nres_base(2,ii)+ipatt(2,i)
5258 energia(kk)=ener(kk,ik)
5260 etot=ener(n_ene_comp+1,i)
5261 rmsnat=ener(n_ene_comp+2,i)
5262 rms=ener(n_ene_comp+3,i)
5263 frac=ener(n_ene_comp+4,i)
5264 frac_nn=ener(n_ene_comp+5,i)
5265 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5266 i,str_nam(ii),ist+1,&
5267 (energia(print_order(kk)),kk=1,nprint_ene),&
5268 etot,rmsnat,frac,frac_nn,rms
5270 write (iout,'(/a/)') 'IEXAM array:'
5271 write (iout,'(i5)') nexcl
5273 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5275 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5276 'Max. time for threading step ',max_time_for_thread,&
5277 'Average time for threading step: ',ave_time_for_thread
5279 end subroutine write_thread_summary
5280 !-----------------------------------------------------------------------------
5281 subroutine write_stat_thread(ithread,ipattern,ist)
5283 use energy_data, only: n_ene_comp
5285 ! implicit real*8 (a-h,o-z)
5286 ! include "DIMENSIONS"
5287 ! include "COMMON.CONTROL"
5288 ! include "COMMON.IOUNITS"
5289 ! include "COMMON.THREAD"
5290 ! include "COMMON.FFIELD"
5291 ! include "COMMON.DBASE"
5292 ! include "COMMON.NAMES"
5293 real(kind=8),dimension(0:n_ene) :: energia
5295 integer :: ithread,ipattern,ist,i
5296 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5298 #if defined(AIX) || defined(PGI)
5299 open(istat,file=statname,position='append')
5301 open(istat,file=statname,access='append')
5304 energia(i)=ener(i,ithread)
5306 etot=ener(n_ene_comp+1,ithread)
5307 rmsnat=ener(n_ene_comp+2,ithread)
5308 rms=ener(n_ene_comp+3,ithread)
5309 frac=ener(n_ene_comp+4,ithread)
5310 frac_nn=ener(n_ene_comp+5,ithread)
5311 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5312 ithread,str_nam(ipattern),ist+1,&
5313 (energia(print_order(i)),i=1,nprint_ene),&
5314 etot,rmsnat,frac,frac_nn,rms
5317 end subroutine write_stat_thread
5318 !-----------------------------------------------------------------------------
5320 !-----------------------------------------------------------------------------
5321 end module io_config